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ABSTRACT 


The object of this project was to define the problem of the assimilation 
of remote sensing data into mathematical models of atmospheric pollutant 
species. An object of remote sensing of the atmosphere is to enable recon- 
struction of the concentration distribution of trace species over a region 
based on the data available from the instrument. The data assimilation prob- 
lem is posed in terms of the matching of spatially integrated species burden 
measurements to the predicted three-dimensional concentration fields from 
atmospheric diffusion models. General conditions have been derived for the. 
reconstructability of atmospheric concentration distributions from data typi- 
cal of remote sensing applications, and a computational algorithm (filter) 
for the processing of remote sensing data has been developed. 
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RECONSTRUCTION OF ATMOSPHERIC POLLUTANT 
CONCENTRATIONS FROM REMOTE SENSING DATA - 
AN APPLICATION OF DISTRIBUTED PARAMETER OBSERVER THEORY^ 


Masato Koda* and John H. Seinfeld 
Department, of Chemical Engineering 
California Institute of Technology 
Pasadena, California 91125 


ABSTRACT 

The reconstruction of a concentration distribution from spatially-aver- 
aged and noise-corrupted data is a central problem in processing atmospheric 
remote sensing data. Distributed parameter observer theory is used to de- 
velop reconstruct! bil i ty conditions for distributed parameter systems having 
measurements typical of those in remote sensing. The relation of the recon- 
structibility condition to the stability of the distributed parameter obser- 
ver is demonstrated. The theory is applied to a variety of remote sensing 
situations, and it is found that those in which concentrations are measured 
as a function of altitude satisfy the conditions of distributed state 
reconstructibil ity. 
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I. INTRODUCTION 

In the remote sensing of tropospheric species, a ground-, aircraft- or 
satellite-based platform performs an instantaneous scan of a region of the 
atmosphere and measures the species burden within the field of view. With air- 
craft or satellite remote sensing the platform is in motion and the field of 
view is constantly changing. An object of remote sensing of the atmosphere 
is to enable reconstruction of the concentration distribution of trace species 
over an entire region based on the data available from the instrument. 

The reconstruction of a concentration distribution from spatially-aver- 
aged and possibly noise-corrupted data is a central problem in processing 
remote sensing data. In the absence of a mathematical model describing the 
spatial and temporal concentration distributions, the reconstruction can be 
carried out by standard data interpolation methods. However, when a mathema- 
tical model exists, the problem becomes one of matching the remote sensing data 
to the model solution in such a way that the incomplete data can be used in 
conjunction with the model to produce an estimate of the region-wide concen- 
tration distribution. This problem of the matching or assimilation of remote 
sensing data into mathematical models for atmospheric constituents is the 
subject of this paper. 

There exist a few recent studies that assess the capabilities of remote 
sensing for monitcring regional air pollution episodes. For example, Barnes 
et al . [1] conducted a comparative analysis of satellite visible channel ima- 
gery in ground-based aerosol measurements. For three cases, each of which 
represented a significant pollution episode based on low surface visibility 
and high sulfate levels, the results show that the extent and transport of 
the haze pattern can be monitored from satellite data. The study demonstra- 
ted the potential of the satellite to monitor both magnitude and aerial extent 


of pollution episodes. In a related study, Lyons, et aU [2] reported on a 
demonstration project showing that currently available synchronous satellite 
data can detect the aerial extent of large scale hazy air masses associated 
with sulfate and ozone episodes. 

A study related to that of the present work was reported by Diamonte, et 
al. [3] in which they considered the comparison of remote and in situ data on 
pollutant concentrations from point sources. They considered typical remote 
sensing geometries to provide insight on estimation of plurne properties from 
these measurements. In a study also related to the present, Kibbler and 
Suttles [4] considered the estimation of unknown parameters in a pollutant 
dispersion model by comparing model predictions with remotely sensed air quality 
data. A ground-based sensor provided relative pollutant concentration measure^ 
ments as a function of space and time. The measured data were compared with 
the dispersion model output through a numerical estimation procedure to yield 
parameter estimates that best fit the data. 

The object of this paper is to define the problem of the assimilation of 
atmospheric remote sensing data into mathematical models of pollutant behavior. 
Since the atmosphere is a three-dimensional system, models of pollutant behavior 
are of the distributed parameter type [5], Remote sensing data represent spa- 
tial averages of concentrations, so that the assimilation problem is, in es- 
sence, one of distributed parameter state estimation. 

First, the concept of distributed state reconstruct bil ity is developed 
for the class of problems of interest. That is, the first question to be 
faced is - can the desired spatial -temporal concentration distribution infor- 
mation be recovered from the measurements in the absence of noise. The deri- 
vation of general conditions that allow one to answer this question is the 
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subject of Section II. In Section III a variety of common remote sensing 
measurement configurations and atmospheric models are tested for reconstructs - 
bility. We conclude in Section IV with general observations concerning the 
inherent potential of remote sensing data in analyzing regional air pollution. 
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II. RECONSTRUCTIBILITY AND OBSERVERS FOR DISTRIBUTED PARAMETER SYSTEMS 

Atmospheric pollutant models consist of partial differential equations, 
linear in the case in which the species does not react chemically or in which 
it is produced or destroyed by a first-order reaction of the form A *►. This 
case represents a wide class of important situations and is the one to which 
we direct our attention here. Nonlinear distributed models must be handled 
by linearization and therefore also fall within the present framework. 

Our interest in this section is to derive distributed parameter observers 
for systems described by linear partial differential equations with inhomo- 
geneous boundary conditions characteristic of atmospheric models. An observer 
is an algorithm that processes measurements of the state of a system to yield 
an estimate of the entire system state. An observer is most frequently 
employed when not all of the states of a system are accessible for measure- 
ment. In the present application, we will be generally interested in only a 
single state variable, the measurements of which have limited spatial resolu- 
tion. The observer is stable if its estimated state converges to the true 
state after a sufficiently long time. The concept of state reconstruct bility 
is useful as a condition for the stability of the observer. Thus, if a meas- 
urement strategy satisfies the condition of state reconstructs lity, then 
the corresponding observer is stable, and, the state (i.e. the concentrations) 
can, in principle, be estimated from the measurements. The condition that 
allows the reconstruction of the system state on the entire field is called 
distributed state reconstruct bility. Associated with distributed state re- 
construct bil ity, the concept of uniform n-mode reconstruct bility can be 


developed. 3oth conditions, n-mode and distributed state reconstructibility, 
will be applied, in Section III, to typical remote sensing measurement 
configurations. 

There exists some previous work on observer theory for distributed param- 
eter systems [6-8]. Kitamura et al . [6] formally extended the lumped param- 
eter observer to the distributed parameter case. Gressang and Lamont [7] 
developed a more complete theory of the distributed parameter observer, includ- 
ing reduced order observers. An application of distributed parameter observer 
theory has been presented by Kbhne [9], The most complete treatment of observer 
theory is that of Dolecki and Russell [8]. In the current work, distributed 
parameter observers are derived in a form appropriate for application to the 
class of systems representing atmospheric species behavior. In addition, a 
result of the present work is an explicit relation between distributed parameter 
reconstructibility and the stability of the observer. Observer stability is 
demonstrated using a technique of Hale [10] in which Lyapunov stability theory 
s extended to function spaces. 

We consider the linear distributed parameter system, 

- L x u(x,t) + B(x,t)/(x,t) (1) 

defined for t > 0, x € D. The domain D is a connected subset of a d-dimensional 
Euclidean space with boundary surface 3D. The d-dimensional spatial coordi- 
nate vector is denoted by x. The state u(x,t) is a scalar function and L x 
is a linear partial differential operator with respect to x. It is assumed 
that the operator L x is well-posed. The input f(x,t) is a known scalar 
function and B(x,t) is a known coefficient. 
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The boundary condition on (1) is 

3 x u(x,t) = h(x,t) X € 3D (2) 

where 3 is a linear, spatial differential operator of suitable order over 
9D and h(x,t) is a known function. The initial condition is assumed to be 
unknown or incompletely known. 

We are interested in considering three types of measurements: 

Case 1: Spatially- Independent Integral Measurements 

The measurement takes the form 

w(t) = / H(x,t)u(x,t)dx (3) 

D 


where H(x,t) is a spatial weighting function. 

Case 2: Spatially-Continuous Measurements 

w(x,t) = C(x,t)u(x,t) (4) 

where C(x,t) is a square-integrable function, i.e., C £ L 2 . 

Case 3: Spatially-Discrete Measurements 

w^(t) = H i (t)u(x i . ,t) i = 1,2,... I (5) 

where w. (t) denotes a measurement at the ith measurement location x^ . By 
taking the limit to small volumes of integration in (3), we can represent 
a system such as (5) by choosing H(x,t) = (t)6(x-x^ ) , i = 1,2,..., £. 
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For the moment let us restrict the problem to one spatial dimension, i.e., 
0 < x < 1. Accordingly, boundary condition (2) can be expressed as 


( 6 ) 


e 0 u(x,t) = h 0 (t) 

3jU(x,t) = h x (t) 

Then the solution of (1) and (6) with initial condition u(x,0) = u Q (x) can be 


x = 0 
x = 1 


expressed in the form* 

A 

u(x,t) - / (r,0;x,t)u 0 (r)dr + | | $* (r,T;x,t)B(r,T)j=‘(r,T)drdT 


:,t) - J 

o 


■IT 

10 


•t A 


ii 

*/ J 
0 0 


$ (r,T;x,t)g(r,T)drdt 


(7) 


where 


g(x,t) = Zb^t)^-!) - 2h 0 (t)6(x) . 


The adjoint Green's function $'(x,t;y,T) is governed by 


(8) 


* ♦ L>*(x,t ; y, T ) = 0 


(9) 


with the terminal condition 

$*(x,t;y,t) = 6(x-y) 

*The explicit form of operators L x , G Q , and ^ are assumed as follows: 

L x (*) = a 2 (x,t) 9 + a^(x,t) + a Q (x,t)(*) 

3X 

6 0 (-) = a 2 (0,t) + e 0 (t)(-) 

Pit-) - a 2 (l,t) + 0j(t)(-) 


( 10 ) 
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and boundary conditions 


* * 

V 


= 0 


( 11 ) 


The operators L x » G Q , and 3j are the adjoints of the operators Bq, and $ 1 
respectively. 

The extension of the adjoint Green's functions to higher spatial dimen- 
sions is straightforward. In higher dimensions, (9) and (10) remain the 
same with the general boundary conditions 


(x,t;y,-r) =0 x £ 3D (12) 

where B is the adjoint of the operator 3 . In general, we note that $ is 

A X 

related to the Green's function $ associated with the system (1) with 
homogeneous boundary conditions by the relationship $(x,t;y,T) = $ (y,T;x,t) 

The adjoint Green's function for well -posed distributed parameter sys- 
tems can be constructed in a variety of ways. Expansion in spatial eigenfunc 
tions and construction of the adjoint Green's function from eigenvalues and 
eigenfunctions is a powerful method for linear systems. Let us assume that 
L has an infinite series of discrete eigenvalues {A-}, i - 1,2,.... Using 

A I 

standard methods, the adjoint Green's function that satisfies (9 )- ( 12 ) is 
found to be [111 

* A -A p (t-T) 

$ (x,t;y,-r) = 2^ <f> n (x)<f> n (y)e (13) 

n=l 

where the eigenfunctions (cp^. } , i = 1,2,..., are the solution of the equation, 
LJ>. = A.<j>. , satisfying the boundary conditions (11) or (12). 

XI II 
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II. 1 Recons tructiblllty Conditions 

The objective of an observer is to reconstruct the system state when 
the measurements are incomplete. To be able to reconstruct the state the 
observer must be asymptotically stable. 

An identity or non-reduced observer for the system (1) with measure- 
ments (4) takes the form 


' 3U 3t^ = L x ^ x>t) + B (*>t)f(x s t) 

+ G[w(x,t) - C(x»t)u(x,t)J (14) 

A 

where u(x,t) is the observer output and G is a suitably chosen integral opera- 
tor with the kernel G(x,y,t). 

Before presenting a derivation of the observer, we will establish the 
conditions under which the system (1) and (4) is reconstructible. We define 
the reconstruct bility kernel function by 

Q(x,y,t) = / / $*(x,t;r,x)C 2 (r,T)$*(y,t;r,T)drdT (15) 

0 u 

It will be shown later that the observer (14) is stable if Q(x,y,t) has a so- 
called generalized inverse, i.e., if there exists P(x,y,t) such that 


/ p 


P(x,r,t)Q(r,y,t)dr = 6(x-y) 


(16) 


By formal differentiation of (15) with respect to time and use of the 
properties of the adjoint Green's function (9) - (12), it is found that 
Q(x,y,t) satisfies the following Lyapunov equation, 
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aQ(x,y t t)_ 

at 


L*Q(x,y,t.) - Q(x,y,t)L* + C 2 (x,t)6(x-y) (17) 


with the initial condition 

Q(x,y,o) = 0 (18) 

and boundary conditions 

f3*Q = 0, Q8* = 0 (19) 

where L*Q = QL*. Although Q(x,y,t) is formally defined by (15), it is impor- 
tant to note that Q(x,y,t) may be computed from (17)-(19) without using the 
adjoint Green's function. 

By using the identity 


8P( 


= -ff P(x,r,t) 3C ^ s,t - P(s,y,t)dsdr 
D J D 


( 20 ) 


P(x,y,t) can be shown to obey the following Riccati equation, 


— 1 L v p ( x >y> t ) + P(x,y,t)L 


■ f P(x,r,t) C 2 (r,t)P(r,y,t)dr 

m/f-K 


(21) 


with boundary conditions 


3 X P = 0, PB = 0 


( 22 ) 


R(x,y,t) may be considered as the kernel of the integral operator P defined as 


The impact of observation error on the design of an observer can be assessed 
from (21) by comparing P to that from the corresponding distributed param- 
eter filter. 
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P/(x) B f P(x,y,t)/(y)dy (23) 

J d 

for f € Lg. 

A linear distributed parameter system (1) and (2) with measurement (4) 
is said to be distributed state reconstmetible if and only if Q(x,y,t) de- 
fined by (15) has a bounded generalized inverse P(x,y,t) for t > 0. It may 
be shown that Q(x,y,t) has a bounded generalized inverse when Q(x,y,t) in 
bounded and positive-definite for t > 0 till . The system (1), (2), and (4) 
will be defined to be uniformly n-mode reconstmetible if there exists posi- 
tive constants M^, Mg, and a such that 


< f/~4> n (x)Q 0 (x,y,t)ct> n (y) dxdy < Mg 


for all t > 0, where <j> (x) is the eigenfunction of L and the modified re- 

n a 

constructibility kernel Q a (x,y,t) is defined by 


x,y,t) = // $*(x,,t;r,T)C^(r,T>^(y,t;r,T) drdx 


t-c D 


The system is distributed state reconstructible if (24) is satisfied for each 
of the eigenfunctions. The uniform n-mode reconstructibil ity test (24) is 
useful when P(x,y,t) cannot be found directly from Q(x,y,t). Since it is 
straightforward to extend the concept of distributed state reconstructibil ity 
to measurement Cases 1 and 3, detailed discussion is omitted here. 

Positive-definiteness of the kernel Q implies that 



f(x)Q(x,y,t)/(y) dxdy > 0 


for all t > 0 and f £ Lg 



- 12 - 


II. 2 Minimum Variance Observers 
11.2*1 Observer for Case 1 

For the system described by (1) and (3), we define the reconstruct!' - 
bility kernel function by 


Q(x,y,t) = f f $*(x,t;r,T)H (r,x)dr f H(s,t)$* (y,t;s,T)dsdr (26) 

4 ) 4 ) 4 ) 


where Q(x,y,t) obeys 


' 3Q ^3t yj ' ^' = " - Q(x,y,t)L* + H (x,t)H(y,t) (27) 


with initial arid boundary conditions given by (18) and (19). Assuming that the 
system is di-otributed state reconstructive, the existence of the general- 
ized inverse P(x,y,t) of Q(x,y,t), that satisfies 


"^ft^ 1 = L v p ( x »y»t) + P(x,y,t)L v 


- I P(x,r,t)H (r,t)dr / H(s,t)P(s,y,t)ds 
T) D 


(28) 


will establish the observer for the system (1), (2) and (3). 

Following Meditch [ 12 ], we define the cost functional associated with 


the observer as 
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■H 


[u(x,0) - u 0 (x)l A x [u(y,0) - u 0 (y)]dx 


j J (w(t) - J H(x,t)u(x,t)dxJ dt 


( 29 ) 


where is an arbitrary final time, u Q (x) is an initial estimate of u(x,0), 
and 


A x (-) 



/ 


P 0 (x,y){-}dy 


-1 


(30) 


P Q (x,y) is a bounded, symmetric, and positive-definite weighting func- 
tion. The observer is found by selecting u(x,t) so as to minimize (29) sub- 
ject to (1) and (2). By minimizing the augmented functional, 


0 = J Q + 


tf 

l 


x(x,t) 


8u(x,t) 

at 


L x u(x,t) - B(x,t)/(x,t) 


dxdt (31) 


the result is the Euler-Lagrange equation, 


3X(x,t) 

at 


L*X(x,t) - 


H(x,t) 



- j H(y,t)u(y,t)dy 
D 


(32) 


with the transversal ity conditions, 

X(x,0) = A x [u(y,0) - u Q (y)] 
X(x,t f ) = 0 


(33) 


Equations (32) and (33) constitute a two-point boundary value problem 
that may be solved by the sweep method. We assume the following Riccati 

A 

transformation for u(x,t), 
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u(x,t) = / P(x,y,t)*(y,t)dy + p(x,t) (34} 

u 

where the kernel P(x,y,t) and p(x,t) have to be determined. 

Substitution of (34) into (1), ( 32 ) and (33) yields 


9p ^ - = L x p(x,t) + B(x,t)/(x,t) 


I 


+ / P(x,r,t)H (r,t)dr 


w(t) 


- / H(s,t)p(s,t)ds 
D 


p(x : 0) = U Q (x) 

6„p(x,t) = h(x,t) x e 3D 


(35) 

(36) 

(37) 


3P(x , Ay . ,t| = L x P(x,y,t) + P(x,y,t)L v 


- J P(x,r,t)H (r,t)drj H(s,t)P(s,y,t)ds 


(38) 


P(x,y,0) = P Q (x,y) ( 39 ) 

6 x P(x,y,t) = 0, P(x,y,t) 3 y =0 x,y £ 3D (40) 

A 

Equations (33) and ( 34 ) imply that p(x,t f ) = u(x,t f ) is the state esti- 
mate at an arbitrary final time t^. It is important to note that (38) is 
identical to (28). Thus we may conclude that the symmetric, positive-definite 
kernel P(x,y,t) completely characterizes the minimum variance observer. 

Equation (35) can be rewritten as 

A 

5 8 t^ = L x “( x > t ) + B(x,t)/(x,t) 

+ K(x,t) I w(t) - / H(y,t)u(y,t)dy 
L D 


( 41 ) 
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where a time-varying observer gain K(x,t) is defined by 


K(x,t) = | P(x,y,t) H(y,t)dy 

0 


(42) 


The structure of the observer is identical to that of the distributed param- 
eter filter [13 ] . 

We introduce the reconstruction or observer error e(x,t) = u(x,t) - u(x,t). 
Then we obtain the following equation for e(x,t), 


— al "—- 8 L e(x,t) - K(x,t) / H(y,t)e(y,t)dy (43) 

L D 

with initial and boundary conditions, e(x,0) = u Q (x) - u(x,0), and 3 x e(x,t) - 0. 
If the initial state is known exactly and the observer is initialized such 

A 

that u(x,0) = u(x,0), then the observer will reconstruct the state exactly. 

It is not reasonable, however, to expect that the initial state will be known 
exactly. It is, therefore, important to insure that if errors are present in 
the initial conditions applied to the observer that the estimate will converge 
to the true value of the state, i.e,, the reconstruction error e(x,t) must 
have the property lim ||e(x,t) j| = 0, for all e(x,0). 

Asymptotic stability of the observer can be demonstrated by using (16), 
(26). (27), and (43). We will consider a Lyapunov function defined by 

V(e,t) = 

It is first necessary to note that Q(x,y,t) is positive-definite and bounded 
from below. Then the time derivative of the Lyapunov function is calculated 
using (27) and (43). The result is 


JJ e(x,^)Q(x,y,t)e(y,t)dydx . (44) 

DO 


i 


u 
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^ V(e,t) = - jj e(x,t) H(x,t)H(y,t)e(y,t)dydx (45) 

DD 

which is a negative-semi definite quadratic form. This is sufficient to show 
that (43) is stable in the sense of Lyapunov [10]. 

II. 2. 2 Observer for Case 2 

In a similar manner to that of Case 1, we can obtain the minimum vari- 
ance observer for Case 2, i.e., for spatially-continuous measurements (4). 

The observer dynamics are described by 

A 

= L x u(x,t) + B(x,t)/(x,t) 


+ / G(x,y,t)[w(y,t) - C(y,t)u(y,t)]dy (46) 

D 

with initial and boundary conditions 

u (x, t) = u Q (x) (47) 

3 v u(x,t) = h(x,t), x € 3D (48) 

A 

where the optimal gain kernel G(x,y,t) is defined by 

G(x,y,t) = P(x,y,t) c(y,t). ^.(49) 


The Riccati equation for P(x,y,t) in (49) is identical to (21) with boundary 

A 

conditions given by (22). The reconstruction error e(x,t) = u(x,t) - u(x,t) 
satisfies 

3e(x,t) 

3t 


= F e(x,t) 


(50) 
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where the integro-differential operator F is defined by 


F e(x,t) = L x e(x,t) 



G(x,y,t)C(y,t)e(y,t)dy 


( 51 ) 


We can demonstrate the stability of the observer by using the reconstruct - 
bility kernel Q(x,y,t) defined by (15) and the Lyapunov function (44). 

Under the reconstructibility assumption, the derivative of the Lyapunov func- 
tion becomes 


which is a negative-semi definite quadratic form. 

II. 2. 3 Observer for Case 3 

For the spati ally-discrete measurements (5), i.e.. Case 3, the observer 
is given by the following system: 

A 

^ = L x u(x,t) + B(x,t)/(x,t) 

z 

+ 22 G.(x,t)[w.(t) - H. (t)u(x. ,t)] (53) 

i=l ITT 1 

where 

G.(x,t) = P(x,x.,t)H i (t) (54) 

and 

• 9P (ft y ' ^ = L x p (*,y.t) + P(x,y,t)L y 
£ 

- 22 P(x,x. ,t)H. (t)H. (t)P(x. ,y,t) 
i=l tit i 


/ 


e (x,t) C^(x,t)e(x,t)dx 


(52) 


( 55 ) 
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with initial and boundary conditions given by (47), (48) f (39), and (40). 

Under the distributed state reconstruct! bi 1 ity assumption, the stability of 
the observer can be demonstrated. 

II. 3 Comments 

The relationship has been established between distributed state recon- 
structibil ity and the existence of an observer. Distributed state reconstruc- 
ts lity is defined through the existence of the generalized inverse to the 
reconstruct! bi 1 i ty kernel. The kernel associated with the observer gain 
satisfies the same Riccati equation as does the generalized inverse of the 
reconstructibility kernel. 

III. REMOTE SENSING MEASUREMENTS AND ATMOSPHERIC MODELS 

In this section we will test both the n-mode and state reconstructibility 
of common remote sensing measurements with models of atmospheric pollutant 
behavior. By far the predominant mode of remote sensing is to measure the inte- 
grated quantity (burden) of material between the ground and some known altitude. 
Thus, both cases we consider here involve vertically integrated data. Various 
assumptions concerning the horizontal characteristics of the measurements will 
be tested. Three-dimensional models of pollutant behavior are generally based 
on the atmospheric diffusion equation [5] that describes the flow and diffu- 
sion of species. The object of this section is to ascertain if the customary 
remote sensing measurements allow one, in principle, to reconstruct the detailed 
concentration distribution. The distributed parameter reconstructibility con- 
dition derived in Section II will therefore be tested in each case. 
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III.l Measurements in a Layer with Horizontal Homogeneity 

The vertical concentration distribution of a pollutant in a layer with 
horizontal homogeneity can be described by the one-dimensional diffusion 
equation, 


3u _ 
3t “ 

subject to 

v iy. 

K 3z 

kIt 


K ^ 


3z 


" h o (t) 

z = 0 

= 0 

z = 1 


(56) 

(57) 

(58) * 


where h Q is a given flux at the ground (z=0) and K is the turbulent diffusion 
coefficient. 

The adjoint Green's function for the system (56) -(58) is 


$*(zj,t;z' ,t) = 1 + 2 X) cos(mrz)cos(mrz' )e 

n=l 


(mr)^K(t-T) 


(59) 


State reconstructibil ity is then to be assessed by condition (24) using the 
modified reconstructibil ity kernel (25). 

We consider each of the measurement types (3), (4) and (5). The condition 
for uniform n-mode reconstructibil ity is (24), which is written for $ = cos(mrz), 
n = 0, 1,2,. . . , as 


-WT 
0 0 


cos(mrz)Q a (z,z' ,t)cos(mrz' )dzdz' < M 2 <°° (60) 


For each of the three types of measurement, the integral in (60) is: 
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Case 1: Spatial ly-Independent Integral Measurements 


/ 

t~CT 


s 2(m0 2 K(t-T)j 


.1 


J H(r,T)cos(mrr)dr| dt 
0 


( 61 ) 


Case 2: Spatiallv-Continuous Measurements 


/' 

t-a 


2(n7r)^K(t-T) 


/ 


C(r,T)cos(mTr)| drdt 


(62) 


Case 3: Spatial ly-Discrete Measurements 


f e 2 - n7r ^ I J (*r) c °s(n7rz^ )} dt 

t-a , ( i=l ’ 


(63) 


for n = 0, 1,2 

From (61)-(63), we see that uniform n-mode reconstruct bility is com- 
pletely dependent on the form of the measurement weighting functions, H(z»t), 

C(z,t) and H.(t) and on the eigenfunction, cos(mrz). The condition (60) 

1 A 

implies that J H(z,t)cos(mrz)dz f 0. We may note that this inequality is 

o 

essentially equivalent to the observability condition derived by McGlothin 
[14]. Similarly, (63) implies that the system state is reconstructible by 
point sensors if the sensors are not located at the zeros of any of the eigen- 
functions. 

In the remote sen c ing problem, the measurement weighting functions are 
often taken as H(z,t) = 1 or C(z,t) = 1. When H(z,t) = 1, the condition (60) 
holds only for n = 0 implying that the spati ally-independent integral 
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measurements do not allow reconstruction of the system state on entire 
•fields.* This can be directly checked by computing the reconstructs lity 
kernel Q(z,z',t). The system with integral measurements cannot be distribu- 
ted state reconstructive since the generalized inverse of Q(z,z',t) = t does 
not exist, since Q is not an explicit function of z and z'. 

When the measurements are spatially-continuous and C(z,t) = 1, the system 
is distributed state reconstructible. From the definitions of Q(x,y,t) and 
P(x,y,t) in (15) and (16), we have 

oo 2 

Q(z,z' ,t) = t + --i- ^ \ cos(mrz)cos(mrz' )(e 2 ^ n7T ^ ^ - l\ (64) 
iTK n=l n L J 


and 


P(z,z',t) = -jr + 4 tt 2 K ^ n 2 cos(mrz)cos(mrz' )|^e 2 ^ nTr ^ Kt - 1 J (65) 

We may note that the integral equation (16) is satisfied when it is recognized 
that 


A mode associated with the eigenfunction 4> =1 (n = 0) can be reconstruct^! 
and the appropriate observer is 

f 1 * 

w(t) - / u(z' ,t) dz' 


A OA 

3u - 1/ + I 

« 3Z 2 * 


Stability of the observer can be demonstrated by constructing the Lyapunov 
function 


V(e 


,t) = J f e(z, 

o J o 


t)te(z‘ ,t) dzdz' 
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00 

5(z-z') = 1 + 2 A cos(mrz)cos(mrz' ) (66) 

n=l 


P(z»z',t) is bounded and positive-definite, and for t > 0, the series (65) 
is uniformly convergent. 


I I I. 2 Measurements of a Steady State, Point Source Plume 

The concentration distribution in a plume from a continuously emitting, 
elevated point source can be described by 


9u _ » 
at " N y 




(67) 


where t is the time an element of fluid spends in the plume from emission, 
equal to downwind distance x divided by the wind speed. The source is of 
strength q located at t * 0, y = 1/2, z * z H (0 < z R < 1), The boundary condi- 
tions on (67) are 


u(0,y,z) = q*5(y - 

1/2)<s(z-z h ) 

(68) 

O 

II 

=J|>» 

CO \CO 

y = o,i 

(69) 

K — <« h 
N z 3z o 

z = 0 

(70) 

|^=0 

3z 

z = 1 

(71) 


The adjoint Green's function for this system is 


C 03 

s>*(y ,z,t ; y' ,z' ,t) =< 1 + 2 53 cos{niry)cos(mry' )e 

V. n=l 


(mr) 2 K (t-t)' 


{ 


» (imr) 'K (t-t) 

1 + 22. cos(m7Tz)cos(rmrz ! )e 


( 72 ) 
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Consider first a scanning measurement performed at a horizontal position 


y B y*» 



J(z)u(t,y*,z)dz 


(73) 


where the scanning data w(t) are taken on a coordinate that moves along the 
t-axis. J(z) in (73) is the altitude-dependent weighting function for the 
measurements. When J = 1, the reconstructibility kernel function becomes 


Q (t,y,zjy',z') = t 


2 y* £Q. s (n^) r cos(mry) + C os(mry' )\fe ^ H - 1 
K H m (mr ) 2 l J V 


cos(n7Ty*)cos(nrny*) 


i-Y.Z , , 

"H n=l m=l (n*+m*)ir 


2^2 v 2 


cos(mry)cos{imry' )< e 


r(nW)ir 2 K H t 


- I 


(74) 


The system is not distributed state reconstructive since the generalized 
inverse of Q(t,y,z;y' , z ' ) does not exist. Therefore, we con' 13 that 
the scanning measurement (73) cannot, in principle, allow rer 't.ion of 
the system state. 

The same results can be obtained for the following measurement systems: 


w(t) = 

w(t,y) 

w(t,z) 


J j u(t,y,z)dydz 


0 0 


= j u(t,y, 


z)dz 


= J u(t,y, 


z)dy 


(75) 


(76) 


(77) 
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In these cases, the reconstruct! bility kernel function Q(t,y,z;y' ,z' ) cannot 
be written explicitly in terms of all the spatial variables y,z, and y',z'. 
Thus, the generalized inverse P(t,y,z*,y' ,z' ) does not exist, which allows 
reconstruction of the system state on the whole domain. As a rule, if 
Q(t,y,z;y' ,z' ) is expressible as an explicit function of all the spatial vari- 
ables and if it satisfies the uniform n-mode reconstructibility test, then the 
system is distributed state reconstructive. 

Indeed, we can show that the system state is distributed state recon- 
struct!' ble for the measurement 


w(t,y,z) = u(t,y,z) 


(78) 


In this case, we have 


Q (t,y,z;y' ,z') » t + K 


v i r t 

> _± — =• cos(nTry)cos(niry' )< e - 1 ( 

H n=l (nit) 2 V. J 


f 2 (n7r) 2 Kyt 

/, — c '(nirzJcosCmTz' ) ( e - If 

n=l (nir; L ** 


00 00 2 2 
+ 2 y y cosCmTylcosttmrzlcosCmT.v 1 IcosCmrrz 1 ) e 2 ^ n71 ^ K H + (m7r) t_ ^ 

n=l m=l (mr) 2 K H +(rmr) 2 K v L_ 


( 79 ) 


The generalized inverse of (79) is given by 
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P(t,y 9 z;y',z') 


1 ^ o f 2(nir) Knt V“1 

~ + 4K h 2 (nffr cos(mry)cos(mry')^e - lj 

S ? /" 2(mr) 2 K v t Vl 

+ 4K V 2 (nir ) cos(mrz)cos(mrz' ) ^e - lj 

8 ]£ ]£ ^(mr) 2 K H + (imr) 2 K v J cos(nTry)cos(rnirz)cos(mty' )cos(nnrz' ) 


2{(mr) 2 K H +(imr) 2 K v }t ~ 


-1 


(80) 


where (80) satisfies the Riccati equation associated with' the measurement (78), 
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IV. CONCLUSION 

This paper has examined the possibility of estimating atmospheric species 
concentration distributions from remote sensing data. Atmospheric concentra- 
tions can be modeled by partial differential equations of the diffusion type. 
Remote sensing data generally represent spatial averages of the concentrations, 
frequently in the vertical direction. The essential problem, therefore, is to 
assess the possibility of estimating the state of a distributed parameter sys- 
tem on the basis of spatially-averaged measurements. The theoretical basis 
of the assessment is a condition for state reconstructibility of distributed 
parameter systems. (The connection between state reconstructibility and the 
stability of the distributed parameter observer has also been developed.) 

A variety of remote sensing measurement configurations were tested for 
reconstructibility. It was found, not unexpectedly, that those measurements 
based on integration of the vertical concentration distribution over the 
entire layer do not lead to distributed state reconstructibility* i.e., there 
does not exist a generalized inverse of the reconstructibility matrix kernel 
and therefore do not afford the possibility of estimating the concentration 
distribution over entire field. Those measurement configurations that, on 
the othe v> hand, enable sampling of the concentration at vertical positions 
lead to distributed state reconstructibility. 
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Alntraet — Optimal filuring itxl trooolhing algorithms for Un«ar dis- 
crete-timr distributed parameter s> stems are derived b> a unified approach 
based on the \N iener-Hopf theon. The Wiener-Hopf equation for the 
estimation problems is derived using the lein-xquare* estimation error 
criterion, Using the basic equation, three types of the optimal smoothing 
estimators are derived, name!;, fixed-point. fixed- interval and fixed-lag 
smoothers. Finally, the refus'd obtained are applied to estimation ol atmo- 
spheric sulfur dioxide concentrations in the Tolushimi prefecture of 
Japan. 

I. Introduction 

A NUMBER of important physical phenomena may be 
modeled as discrete-time distributed parameter sys* 
terns. When estimation problems are encountered in such 
systems, the measurements are also frequently discrete in 
time. A great deal of work has been carried out on estima- 
tion problem* for continuous-time distributed parameter 
systems P]-|4). Tzafesta* [5). [6] and Nagaminc et a!, [7] 
have derived optimal estimators for discrete-time distrib- 
uted parameter systems. Tzafcstas employed a Bayesian 
approach, where Nagaminc et al. considered only the filter- 
ing problem based on the Wiener-Hopf theory . Recently 
Bencala and Seinfeld [3] have derived the optimal filter for 
continuous-time distributed parame«er systems^ with dis- 
crete-time observations by the Wiener-Hopf approach. 

The object of this paper is twofold. First, we seek to 
derive optimal filtering and smoothing algorithms for dis- 
crete-time distributed parameter systems by unified 
Wiener-Hopf approach. Fixed-point, fixed interval, and 
fixed lag smoothers are considered Second, we wish to 
apply the results to the estimation of atmospheric sulfur 
dioxide concentrations in the Tokushima prefecture of 
Japan. 


II. Description of the Distributed Parameter 
System f 

Let D be a bounded open domain of an r-dimensional 
Euclidean space with smooth boundary 9 D, The spatial 
coordinate vector will be denoted by x ~ (x j,* * v, x f ) 6 D, 
Consider a linear distributed parameter system described 
by 

u(/c + 1, x) — E x u(k, x) + 6(Ar, *)*•(&, x), x & D 

0 ) 

where u{k 4- l,x) is an n-dimensional vector function of 
the system, h*(A\ x) is a vector- valued Gaussian process, 
is a linear spatial matrix differential operator, and G(k, x) 
is a known matrix function. 

The initial and boundary conditions arc given by 

u(0. x) = u 0 (x) (2) 

r { u(*+l. {) = $(*+ 1,0. iedD (3) 

m-] = «t-]+(l-«(0)3l-)/3" (<) 

where n is an exterior normal vector to the boundary HD at 
a point f £ 3 D and o(0 is a function of class r 3 on 3 D 
satisfying 0 < a{£) < 1. S(k + 1. f) denotes a source func- 
tion at the boundary' and is assumed to be known 
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Assume that u$(.v) is a Gaussian random function the 
mean and covariance functions of which are given b> 

£ K(*)] = ° ( 5 ) 

£[« 0 (*)woO')) = Fo(*' >') (6) 

where £[*) and the prime symbol denote the expectation 
and transpose operators, respectively. 

Let _ihe observed data be taken at m points. ***♦*•♦ 
x m 6 P *= D U 3 D and let an mu-dimensional column vec- 
tor ujk) be defined b> 

u m (k)=-Co\[u(k.x')."-,u(k.x m )}- ( 7 ) 

Let th£ observations be related to the states by 

:(k)~H(k)u m (k) + v(k). (8) 

where :(k) is a p-dimensional observations vector at the m 
observation^ points, xK* • % x m e D. H(k) is a known 

pxmn matrix, and v(k) is a p-dimensional vector-valued 
white Gaussian process. Assume that the white Gaussian 
process u*(A%a*) in (11 and v{k) in (8) are statistically 
independent of each other and also independent of the 
initial condition u 0 (x). Their mean and covariance func- 
tions arc given by 

£[h(A\*)] = 0. £[c(A)] =0 (9) 

£[«•(*:. jt)M’(j.v)] = Q(k. x, y)S it . x. y 6 D 

(10) 

E[v(k)V(s)} = R{k)S kl , (11) 

where B kl is the Kronecker delta function, and Q(k . x. y) 
and R(k) are s>mmctric positive-semidefinite and positive- 
definite matrices, respective!). 

III. Description of the Estimation Problems 

The general problem considered here is to find an esti- 
mate u( t. x, k ) of the state u ( ?, x ) at time t based on the 
measurement data r£. denoting a family of r(o) from 
o = 0 up to the present lime k. Specifically, for t > k we j 
have the prediction problem, for t = k the filtering prob- j 
1cm, and for r <k the smoothing problem. As in the 
Kalman-Bucv approach, an estimate u(r. x/k) of u(t, x) 
is sought through a linear operation on the past and 
present observation values zjj as follows! 

k 

u(t.x/A)= 2 £(T.,\\cr)r(o) (12) 

o = 0 

where /(r, x o) is an nxp matrix kernel function. 

To differentiate between the prediction, filtering, and 
smoothing problems, we replace (12) with different nota- 


tion for each problem: 

1) Prediction (r > k) ( 

u(t,x/*)= 2 ^(T,x t c)e(o). (13) 

2) Filtering (r = k) 

d(k,x/k) = 2 F{r,x,a):(o). (14) 

0 = 0 

3) Smoothing {r < k) 

k 

,x/k)~ 2 £( T ,*.Xc)x(o). (15) 

o = 0 

The estimation error is denoted by u(t, x/k), 

6[r,x/k) = u(7,x)-iHT.x/y). (16) 

The estimate fi(T. x/k) that minimizes 

J(u) ~ £[lifi(r. x/A)l |3 ] (17) 


is said to be optimal, where 1! * 1! denotes the Euclidian 
norm. 
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Theorem i (Wiencr-Hopf theory) A necessary and 
sufficient condition for the estimate i3(t, a f k) to be opii* 
ma! is that the followjng^Wiener-Hopf equation holds for 
a =• 0, !*••*. k and a* G D . 

2 At. x. o)£[;(c).-'(a)] = £[u(t, x).*'(a)]. 

»£>C 

(18) 

Furthermore. (IS) is equivalent to 

£[i?{t, x/A );'(a)] = 0 (19) 

for a = 0, L« • », Ar and A G />. 

/V 00 / Let be an tixp matrix function and 

let c be a scalar-valued parameter. The trace of the covari- 
ance of the estimate* 

tf«(T,x/A')= 2 (F(T*x.o) + <F A (7.x t o)):{o) 

o=0 


is given by 

A*,) = e\ 


hu(r.r)- u(t, x/k) - * 2 ^(t. a:. o)r(o )** : 

OB(* 

k 

u'(i.x/k) 2 /i(r.x*o)r(o) 

o = 0 


= £[r fi( t. *.•*)! : ] -2tE\ 

k 

+ r 2 £[ 


2 ^(T*A\o):(o)f: J j. 

A necessary and sufficient condition for u(t. x/k) to be 
optimal is that 


dt 


= 0 , 


that is. 


w'(t, x'k) 2 f A (T. A*o)r(o) 


= 0 


for any nxp matrix F^ir. a*, a). Using the relation between 
the trace and inner product yields 


B'(r.xA) 2 FJr.x . o).-(o) 


= tr 


fi( T.x/k) 2 S'(o)A( T.AT.o) 


= 2 tr [ £[ u( r. x /A, )r '( o )] t. X. o )] = 0. 

0 = 6 


Setting a* A) = £[w(r, a /A) r'(o)) in the above equa- 
tion, it follows that (19) is a necessary’ condition for 
u(r.x/k) to be optimal. Sufficiency of (19) also follows 
from the above equation. Q.E.D. 

Corollary l: (Orthogonal projection lemma). The follow- 
ing orthogonality condition holds, 

£[fi( T ,x/A-)fi'(f. .>•/*)] = 0, x, y € D (20) 

where £ is any time instant, for example, £ < k. £ — k or 
i>k. 
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Proof Multiplying each side or (19) bs As* X * <*) and 
summing from o s= 0 to a ~ A' yields 

k 

£ fi{r.*'A) 2 r'(a)£'(>\,v.£t) = 0, 

•«o 

Substituting (12) into the above equation vields (20). 

Q.E.D. 

Then the following lemma can be proved. 

_ Lemma l (Uniqueness of the optimal kernel). Let 
At, *,o) and At, x, o) + A'(t, x, p) be optimal matrix 
kerne] functions satisfying the Wiener- Hopf equation (18). 
Then it follows that 

A’(t, jc, o) e 0, <r = 0J , * • ♦ , A and x € D> (21) 

and the optimal matrix kernel function A t, a. o ) is unique. 
Proof From (18) we have 

2 F(r,x, o)£[r(o)r'(o)] 

6 5=0 

= £[«(T.A)r'(o)] 

= 2 (F(r.x,a) + A'(T.A.o))£[(r(o);’(a)]. 

o=0 

Thus. 

k 

2 A’(r. a, o)£[r(o)r'(a)] =0. 

0 = 0 

Multiplying each side of the above equation by A''(t. a, q) 
and summing from a = 0 to a = k yields 

k k 

2 2 A’(t. a*. o)£[r(o)r'(a)]A*'(T* x< a) — 0. 

0=0 6=0 

On the other hand, from (8) and ( 1 1 ) we have 

£[ .-I o ).-'(<»)] = Hlo )£[« Jo K,(o)] H‘la) + F(o)S co . 

Then it follows that 

/. k 

2 2 K{?*x.e)fi{o)E[uJo)u , m {a)]H r {a)S , {7.x,a) 

o=0 o-V 

k 

+ 2 F(*.x.o)H(o)R(o)H , {o )A”( t. a. o ) = 0. 

o = 0 

Since both terms on the right side of the above equation 
are positive-semidefinite because of the positive-definiteness 
of /?{o). a necessary and sufficient condition for the above 
equation to hold is A*(t. a, o) h 0. o = 0. !,* * *, k and x & 

D. Thus, the proof of the lemma is complete. Q.E.D, 

In order to facilitate the derivation of the optimal esti- 
mators. we rewrite (18) in terms of the following corollary . 

Corollary 2; The Wiener-Hopf equation (18) is rewritten 
for the prediction, filtering, and smoothing problems as 
follows. 

1) Prediction (r > k) 

2 .4(t, a, a )£[;(o );'(<»)] = £[i/(t, a):'(<x)]. 

6 = 0 

(22) 

for a = 0, 1,' • *, k and x 6 £, 

2) Filtering (t = k) 

2 £(Ar,A,o)£[r(o)r'(a)] = £[u(A, A)r'(a)] (23) 

o=0 

for a = 0. 1.* - *, k and x E D. 

3) Smoothing (r < k) 

2 5(t, k, a, o)£[r(a)r'(a)] = £[u(t. xV'(a)] 

o = 0 ' 

(24) 

for a - 0. k and x € D, 

In what follows, let us denote the estimation error covari- 
ance matrix function by P( r. x, y/k), 

P(r,x, y/k) = £[c(t. a/A‘)u'(t. y/k)}. (25) 
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III Derivation of the Optimal Predictor 

Jn this section, we derive the optimal prediction estima- 
tor b> using the Wiener-Hopf theory in the previous 
section. 

Theorem 2: The optimal prediction estimator is given by 

&(k+\.x/k) = £Ak**/l<) (26) 

16 00. (27) 

Proof; From (22) and (1) we have 
2 /f{A4- l.x.o)£[r(o)e # (a)] = £ s E[ulk. x ).-'(«)] 

c«0 

since h ik. x) is independent of :(a), a ~ 0, 1,**% A\ From 
the Wiener- Hopf equation (23) for the optimal filtering 
problem we have 

2 {/f(A + 1.x, o) - £ x F(A.x,o))£[e(o).-'(a)] = 0. 

o»0 

Defining N{k\ x, o) by 

A' (A + 1. x,a) = A(k + I, jr.jf) — L x F(k>x,o) t 

it is clear that A{k + 1* x,o) **- N(k + I, x, a) also satis- 
fies the Wicner-Hopf equation (22), From the uniqueness 
of A(k ’t U.o) by Lemma 1 it follows that S{k + 
U,o)50, that is, 

A{k+ 1. x, o) = £ x F{k\x. A). (28) 

Thus, from (13) and (14) we have 

k 

u{k + l. x/k ) = L, 2 £(o, x, o)r(o) = l x u(k, x/k), 

okO 

Since the forms of T ( and St A * L€) arc known, the 
predicted estimate t/(A -*-!.{ A) also satisfies the same 
boundary condition as (3), T^uik •+* 1. f /k) - S(k + L {)« 
( € dD Thus, the proof of the theorem is complete, Q.E.D. 

Theorem 3; The optimal prediction error covariance ma- 
trix function P{k t \,x. y/k ) is given by 

P{ A -r 1 . X , V /k ) = £* P( A . x. y/k )£; + Q(k, x. y). 

(29) 

r t P(k+ \.i.y/k) = 0. *63 D (30) 

where 

<?(A. X. y) = C(A. x)e(A.x, y)C‘(k. ,»•). (31) 

/Voo/ From (1). (16). and (26) it follows that 
fi(A + 1. x/A) = U x fi(A. x/A) + G(A. x)vt'(A\ x) 

(32) 

and from (3), (16), and (27), 

r,fi(A+ l,*/A)-0, *6 32). (33) 

Then we have from (31) P(A + 1. x, y/k) = £|u(A + 
1, x/k)u'(k + ). „v/A )) = t x P(k, x, >-/A)£; + jJ(A, x, ,v) 
and from (33). £[r { fi(A + 1, */A)C'(A + 1, >•/*)] = 
T t P(A -f 1, {. y/A ) = 0, Thus, the proof of the theorem is 
complete, Q.E.D, 
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IV, Derivation of the Optimal Filter 

Let us derive the optimal filter by using &he Wiener- Hopf 
theorem for the filtering problem. From (23) it follows that 

F(k + 1, x. A + l)£[c(A + l)r'(a)] 

+ l,x,o)£[r(a)c'(a)] 

0*0 

= £[u(A + 1, x)e'(a)] 
for a = 0, !.•••, A -r 1. 


( 34 ) 


UMlttl** 


From U) and the independence of and *ik * LaK 
it follows that E\u(k + \. x):U)]~ £,E{u(k. x):'{a)}. 
Applying the Wiener- Hopf equation (23) to the right side 
of the above equation yields 

E[u(k+ \.x):-(a)] = £, 2 m.*,°)£[.*( 0 )r'(a)]. 

e*0 

(35) 

Furthermore, from (8) and the whiteness of v(k ■«+ 1) we 
have 


E[:(k + !).•'(«)] - H{k + 1 )E[u m (k + 1 );'(«)]. 
Let us introduce £4*1 ? n< * [*J£* as follows. 

[M-l o 


e.N- 


and 


0 £*-[•] 
I-))'- 


(36) 


(37) 


Then from (1) and (7) it follows that 

u m {k- \) = £.ujk)-r* m {k) (38) 

wjk) - Co\[c(k, x')w(k.x'),---, 

$G{k,x m )w(k.x”')]. (39) 

Then we have for o < A + 1, E[:{k + l):'(a)) = //(/; + 
?}£.£( u ttI ( A )r*< a)} Applying the Wiener- Hop f equation 
(23) to the right side of the above equation yields 

£[£(£+ !):'(•)] 

= tf(*4 l)u 2 •••<’ )£[-•(« )-*'(«)] 

e = 0 

(40) 


where 

Fik.x'.o) 
F(k,x m , o) 
Subsmunng (35) and (40) into (34) yields 


FJk. o) = 


(41) 


J AiU.x.o)£[r(o);'(a)] = 0, o = 0, ),•••, k 
o=0 

where 

AiU.x.o ) = F(/+ 1.x. * + 1)W(*4 1)uF m (*.o) 

a.o) + £(/; + L a.o). 


Since it is clear that F(A\ a. o) + A^(A\ a. o) also satisfies 
the Wiener- Hopf equation (23), it follows from Lemma 1 
that A*j( A\ a. o)5 0, Thus, we have the following lemma. 

Lemma The optimal matrix kernel function F{k< a, a) 
of the filter is given by 

F(k+ L a, a) — t x F(k, a, o) 

-F{k + L a, A* -r 1 )//(£ + \)£ 9 F m (k*o). 
0- o.l. •••,*. (42) 


Theorem 4; The optimal filtering estimate u(k,x/k) is 
given by 

u(k 4 1, x/k 4 1) = £j>(k t x/k) 

+ F{k 4 1.x. *4 l>(*4 1) 

(43) 

»(*■ 4 1) = *(k 4 1) - H(k 4 \)ZJ m (k/k) 

(44) 

u(0. x/0) = 0 (45) 

r t u(A 4 1 . t/k 4 1) = S(k 4 I. f ), { € 3Z> (46) 

where 


u w «./*) = Col[i(A-,x'/A).-... 
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ju 


X.c.k 


*(k.x m /k)]. (47) 
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Proof, Using < 14) and (42) yields 
u(k + hx/k + 1) » F(k + It x, U + 1 ):(k 4* )) 
k 

+£, 2 F(k.x,o)i(o) 

«»0 

-£(* + l,x.* + I )//(* + 1)£. 

• 2 FJk.p)z(o). 

B«0 

Again from (14) we have 

I^/A'4- !) = £,*(*•*/*) 

+ F(k+ Lx.h* ! )v[k+ 1). 

Since wc have no information at the initial time, it is 
suitable :o assume an initial value of u{k 4* 1, x/k 4* 1) as 
w(0, x/0) - £( i/ 0 (at) 1 » 0, Furthermore, since we know the 
exact forms of T ( and 5(A*4- Uf), the boundary value 
u{k + 1. {/A 4- 1) also satisfies the same boundary condi* 
tion as u(k •+• 1. £). Thus, we have T ( u{k + U {/A + |) = 
S{k + l,f)« {£?£* And the proof the of the theorem is 
complete. ~ Q.E.D. 

Note that HA + 1) defined by (44) is rewritten by using 
the prediction value of (26) as follows. 

*(A-r 1) ~ r(A 4- 1) - H(k 4- l)u w (A 4* 1/A) (48) 


or 

p(k 4- 1) as //(A *r !)<?„(* 4- 1/A) + e(A 4- l) (49) 
where 

OJk * 1/A.) = Co)[i5(A + 1. x'/k),--’, 

«* + l.y/*)] (so) 

and 

SJk~ 1 A ) = MA- D-i/JA-r 1/A). (51) 

HA -r 1) is termed the innovation process [8j. [9], 

In order to find the optimal matrh kernel function 
F(k J.,\, A -*■ J) for the filtering problem, we introduce 
the following notation. 

pj t. x 'A ) « [ p( r* x, x'/Ar)** * % p( t, x, x m /A )] 

(52) 

and 


PmJ 7 /k ) 


p(r, x\'k ) 

p( 7. x m /A ) 

f p(T.xKx'/k)r">p(T,x\x m /k) 


[p{T t X m ,X } /k) t "'sP(T,X m .X m /k) j 

(53) 

Note from the definitions of p m (r % x/k) and p mm (r/k) that 
P m (r> x/k) = £[u(t, x/k)V m {v/k)] (54) 

and 

p„Jr/k) = E[a m (r/k)a‘ m Wk)]. (55) 

Furthermore, we define the covariance matrix of the in- 
novation process HA 4* 1) by T(A 4- 1/A), 

T(A + 1/At) = £[r(* + 1 )f'(k + 1)]. (56) 

Then from (49) it follows that 

r(*+ \/k) = H(k+ \)p„„(k+ l/k) 

•H\k + !) + *(*+ 1). 

Then the following theorem holds, 


£ 
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Theorem f The optima) filtering gam mam* function 
Fik + I. a *. k •+ Mu gi\cn b> 

F(k + ).*» A + i) b p m ( A * h^Jt) 

*HU+\)T m Hk + m)m 


or 

F(A-r U* # A + 1)=/UA- 1,*/A) 

**(* + l/A)//'(A + 1)/? W, (A* + 1) 
(59) 

where 

Hk + I A) = (/ + A{k + I )/>„„(* + l/k))~' (60) 

and 


A(k + I) b HU + l )/?“'( A’ + 1)W(A’ + 1). (61) 

Froo/. From the Wiener- Hopf equation (23) it fol- 
lows that 

F(A- + l, *. A + 1 )£[i(A + )):*(*+ 1)] 

k 

+ 2 £(* + Ux.6)E[:(o):U+ 1)] 

9P0 " 

= £[h(A* 1, x):'(A * 0), 

Substituting (42) into the above equation yields 


F(k + 1,*,* + l)|£ 
A.o)r(o )]:'( 


:(A + !)-//{*+ l)t. 


* 2 /*,{A>o)r(o )Jr'(Ar * 1)| 

= £ ■ u{k + 1.x) - 2 £(*,*, o)r(o)} 

1 1 «=o 

Substituting (14) into the right side of the above equation 
and using (26) and the orthogonality condition of (20) 
yields 


)}-•'(* 


+ >) 


\u(k- \.x)-i t 2 f(A\*.o)r(o )}.•'(* + D 

l **< \ 

= £[u(A * 1 ,x/k\;U+ D] 

= £[C( A 1. x/k)u‘ m (k ■+• l)]// - !, < + 1) 

- p fl (k tI,j 'A)H'(A ■+ 1). 


Using the orthogonality condition of (20) gives 

£[KA- + l ).-*(* -r ))] «mk + \)E\ajk + i/lr)fiN 
• / (k* \)]HU + 1) + ~A(k + 1) 

= HU * ))/>„,., (A + 1/Ar) 

•A’(* + 1) + R{k + 1) 

= HA + \/k). (62) 


Then we have 

F(A + \.x,k+ l)r(A+ 1 /k) 

= p m {k+\.x/k)H%k+\). 

(63) 

Thus (58) is derived. In order to show the equivalence 
between (58) and (59), we use the following matrix in* 
version lemma, 

PH\HPH ■ - R)"' = P{1 + H'R-'HP)~'H‘R-'. 

(64) 

Frorj (58) and (64) we have 

£(* + !.*,* + \) = p m (h + \,x/k) 

*+(A j + \/k)H'(k+l)R~'(k+ 1). 

Then (59) is derived, and the proof of the theorem is 
complete. Q,E.D. 

The equation for the optimal filtering error covariance 
matrix function p(A' + 1, x, y/k -r* 1) now must be de- 
rived. 
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Theorem 6 The optimal filtering error covariance matrix 
function p(A t- I. x. t'A + I) is given by 
plk + l.x.y/A + I) >=p( A * 1, a, y/k) 

-/>„(A+I,x'A)/m + I) 

• r“'(*+ l/A)//(A+ 1) 
•/>;<*+ hy/k) (65) 


or 

p(k + l.x.y/A + J )-p(k * l,x,y/A) 

-P m (k+Ux/k)Hk~\/k) 
■Mk +\)p' m ik + l.y/A) 

( 66 ) 

where 

"(Q. x, y/0) =p 0 (*.y) (67) 

and 

r t p(k ■*> l,{, y/A + 1)50, ( e 9/>. (68) 

Proof From (!) and (43) we have 
fi(A + l,x/A* 1) *• t!(A 1.x/ A) 

-£(*•*• 1.X.A + l)K*+ I) 

( 69 ) 

and from (3) and (46) 

T t B(k 1. ( /A + I) = 0. (e 3J>. (70) 

Using the independence property between t’(A + 1) and 
u(A * 1. x/k) or fi(A + 1, y/k ) yields from (69). 

p(k + l.x, y/k + !) = £[«( A * 1.XA + I) 

•u’(A + l,.i/A -r ))] 

~p(k~ l.x.y/A) 

+F(k- I. .v, A -r 1 )£[•>( A ■+• 1) 
•i*'(A * l)j£!A + l.y.A-t 1) 
-£(A + 1.x. A - 1)//(A + 1) 
•£[fi„(A* 1 'A)fi'(A- 1 , y/k )] 
-£[fi(A<+ \,x/kW m (k + 1/A)] 
■//•(Af \)F(k + 1. y, k + 1). 


Using (58) and (63) it follows that 
p(k * I. x, y/k + l)=p(A l,x,.v/A*) 

~Pm(k t- \,x/k)H'(k + I) 

■F(k-r 1, .v. A + I) 

= P(A + I, x, y/k) 

-pj A-M,x/A)H'(A+ 1)0 
.(a + I/A)//(A + l)p;(A +1. y/A). 


Thus (65) is derived, The equivalence between (65) and 
(66) is easiK shown by using (64), Since the initial value 
t/(0, x/0) of u(k * 1, x/k + 1) is zero from (45), it is clear 
that /( 0, jt, y/0) = £[5(0. y/ 0)) “ fax, y). Multiplying 
each side of (70) by u\k 4* 1, y/k + 1) and taking the 
expectation yields T^flk + 1, t y/k + 1) = 0, ( e 9Z>, 
Thus, the proof of the theorem is complete. Q.E.D. 

Corolla o 3 . u m (k + \/k + 1) and + hx/k+ 1) 
satisfy the following relations, 

ojk + 1/A + 1) 

= u„(A + 1/A) + Fjk + 1. A + l)i'(A + I) 

(71) 

F m (k+ l.A+1) 

= Pmm(k+ 1/A)>KA+ 1/A)//'(A+ 1)£'\A + 1) 

( 72 ) 


or 

£„( A + 1, A + 1) = p„„(A -r 1/A + 1) 

•//'(A + l)£" l (A-r I) (73) 
Pmm(k + 1/A + 1) = p mm (k + I/A ) - p mm (k + 1/A ) 

•*(A + 1/A)£(A + \)p„Jk + 1/A) (74) 
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or 

p„Jk - hk 4 I )o Pmm lk + \ 'km 4 1 /kb 

(75) 

/V 00 / From the definitions (41) and (50) of F m (k 4 
I, A *r 1) and <t m (k + 1/A)* it is clear that (71)* (72), and 
(74) hold, From (60) and (74) » follows that 
P mm (k 4 \/k + I)b;J* + 1AW(* + l A) 

* 4- !/A*) - if (A 4 1) 

•Pmmik + *A)J 

«/U* + !/*)«* + »/*) 

.{/ + au + l)A,Jfc+IA) 

* l)Pmm(h + I A)) 

Thus. (75) is derived and (73) is clear from (72) and (75), 

Q.RD. 

The present result corresponds to that of Santis tt ai 
1 17) which is an abstract form of the filter. 

V. Derj\ atios or the Equations for the 
Optimal Smoothing Estimator 

In this section* we derive the basic equations for the 
optimal smoothing estimator b> using the Wiener-Hopf 
theory* 

Lemma 3 The optimal matrix Kernel function B{t, k 4 
1* x* a ) of the smoothing estimator is given by 

B[ r, A -*• Ux.o)*B(uk.x.c) 

-3(T,Hl, jf»H \)H{k 4 l)£.f m (A\d), 
o«0. A. (76) 

Proof From the Wiener-Hopf equation (24) we have 

2i(r,Jl+l,*,«)£[;((i):'(o)]=£[#(T,*|:'(«)], 

0*t' 

o = 0,* • *, A 4 ) 

and ! 

* ! 

2 £(t. A, A\o)£[:(o):Xa)] « £[u(r. x)r'(a)]* 

o**0 

a = 0,* * *, A. 

Subtracting the latter equation from the former yields 
£{?*A 1.x* A 4 1)£[;(A * 1):'{a)] 

A 

4 2 (B(t, A 4 l.X.ff) 

— fi(r, A, x.o))£[;(o)e'(a)] =0* 

From (8) and (23) we have 

B[:{k + ))-”(<•)] = H(k + l)£,£[i/„(A)x'(o)] ' 

= H{k+\)£. 1 FJk.o) 

®*0 

•£[*(o)x'(«)]. 

Then it follows that 
k 

2 A'(t. k, x, o )£[-•(» ).*'(“)] = 0 
®«o 

where 

tf(r, A, x, o) = £(r, A 4 1, x, a) “ B(r, A, x. o) 

4£(t,A 4 1. x, A 4 1)//(A 4 l)t # F m (A, o). 

Since it is easily seen that £(r. A* x, o) 4 $( t* A* x* o) 
also satisfies the Wiener- Hopf equation (24). from Lemma 
1 we have X f (r, A, x, a) s 0* and the proof of the lemma is 
complete. Q.RD. 

Theorem 7; The optimal smoothing estimate u(r, x/k 4 

1 ) is given by | 

w(t» x/k 4 1) = u(t. x/k) 4 £(r, A 4 1.x, A 4 1) 

•'(* + !) (77) 

r c fi(r.t/A*4 l) = 5(r.4).( €0D (78) 

A~ T,T + 1.*'** 
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Proof From 05) it follows that 
w(t, a k ** 1) *= B(r t k + I* jc, k 4 )):{k 4* 1) 

i 

+ 2 B{r.k+ I.Jt.oktV), 

***0 

Substituting (76) into the above equation yields 
&{t,x/k 4 )) 

~£( T .A4 Ujc.it 4 I) 

* ji(A + l)-//(A + I)£. 2 F.(*.*)*(o)J 

A 

4 2 J3(t» A. x*© ):(<?). 

«»o 

Substituting (14) and (15) into the above equation yields 
«(t.x/A4 

4£(r.A4 Ux.Ai* ))KA 4 1), 

Since we have no additional information about the 
boundary value of u(t. a). except for S[ t, £) and the exact 
form T ( * we have £/A 4 1) » S(t, £), £ € 3/3. and 
the proof of the theorem is complete, Q,E D. 

Theortm 8 . The optimal smoothing gain matrix function 
B{ t, A* 4 1, a*. A 4 1) is given by 

B( T. A + I,*, A + 1) - £„( T, x/A)£^£) 

>U + ljr-'(A+ I/A) (79) 
or 

£(r,A-i- l,x.A + 1) = ./(t.x/A + 1) 

•W'(A+ l)«-'(A+ 1) (80) 

where 

J(r. x >k ■+■ I) = L m (r, x/k)£i 

• (/ + fi{k * \)p m Jk - 1 'A)) 1 

( 8 !) 

I m (T. x,'A) = [£(t.x. x'/A ).'••. Mr.r,*"/*)] 

(82) 

and 

i(r. x. >/A) = £[u(t. x/A)C'(A, >-/A)]. (83) 

Froo/; From the Wiener- Hopf equation (24) st fel- 
lows that 

a(r,H 1* a, A* 4 l)£[r(A4 1);'(A*4 1)] 

A 

4 2 ^( T » A 4 1. a. o)E[z(o): , (k 4 1)] 

= £[ U (r,x):'(A+l)]. 

Substituting (76) into the above equation yields 
2 ?{t, A + l,x, A + l)£[v(A + 1).*'(A+ 1)] 

- r[fl( T . x/A).-'(A + 1)]. 
On the other hand, from (4P * and (49) we have 
£[»(A + 1);'(A + I)] = £^(A + 1)(»>(A + 1) 

+«( .A + 1JMA +!/*))'] 

— E[r(k + l)v'(A + 1)] 

= T(A + 1/A). (84) 

From (8) and the independence of c(A 4 1) and 
u( t, a/ A* ), we have 

£[c(t, x/A)»'(A + 1)] 

- £[fi(r, x/k)u' m (k + l/k)]H'(k + 1). 
But from (26) and (38) it follows that 

*Jk + 1/A) = L.0„(k/k) + * m (k). (85) 

Then we have 

B(t, A + l.x.A+ 1)T(A+ 1/A) 

= £„( r.x/A)£iff'(A+ 1). 

1 

4*. 
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using, t he maim inversion lemma (64 1 Thus the proof of 
the theorem is complete, Q.E.D, 

Let us now denve the equation for L(r,x t y/k + 1), 
Using the orthogonahts condition (20) yields 

Z.(t, x. y/k + I) = £[u(r, x)B'(k + \.y/k + 1)]. 
Substituting (69) into the above equation yields 
Hr.x. y/k+\) = L{r. x>y/k)i; 

-L m (r.x/k)\Z l .H'(k + \)F{k+ + 1). (86) 

From (3) and (78) it follows that r ( u(r. t/k + 1) ® 0, f € 
0/3, Multiplying each side by w'(A + l, y/k + 1) and tak- 
ing the expecution yields 

r t t(i\{.,v/*+ i) = o, fea/3. (87) 

Then the following theorem holds. 

Theorem 0 J(l x f k & 1) in (80) is given by 

J ( T. x/k *+ 1 ) = J { r , xjc )£ W (A + 1/A*) (88) 

y(T.x/T) ~pjr.x/r) (89) 

r t /(T.{/* + i) = o. «ear». (90) 

Proof From (86) and (59) it follows that 
k„(v. x/k + )) = L„(t,x)V. 

• {/ - J?(* + + \/k) Pmn (k + 1/A')). 

But we have 

1- R{I~PR)~'P = RU + PR)~ l 

■ ((I + PR)R~' - P) 

= ((/ + pr)r-')~'r~' 

= (7 + APf '. 

Thus. 

LJr.x k*- l) = JJr.x)C: 

■{J-R{ k+\)p m Jk+\/k))~'. 
Therefore, from (8)) it follows that 

At. x/k- 1) - A m (r, x/k + 1) (91) 

and from (81) we have 
J{t. x/k + \) = At , x/k)£‘. 

•(/ + *(* + l)p mm (k+l/k))~\ 

Then it follows that 

At. x/t) = J.Jt.x/t) = pjr, x/r). 

Since (90) is clear from (87) and (91)* the proof of the 
theorem is complete, Q.ED, 

Let us now derive the equation for the optimal smooth- 
ing error covariance matrix function p(r. x, y/k) defined 
by 

p{r , x, y/k) = £[u(r f z/k)u(r> y/k)]< (92) 

From (77) and (78) it follows that 
u(t* x/k + 1) = fi(r, x/k) 

+ 1*x, * + \)p(k + I) 

(93) 

r t fi(T.«/fc+ 1) = 0, {€0/3, (94) 

Then the following theorem holds, 

Theorem JO: The optimal smoothing error covariance 
matrix function p(t, jc, y/k + 1) is given by 

p(r.x.y/k + l)=p{T,x,y/k) 

■{k+l/k)H(k+ l)£ir;(r. y/k ) 

(95) 
or 

pi T. X. ,v/A + 1) = p(t. *. >'/*) 

- J(r.x/k + 1)(0 ) 

• *(* + 1 /k )R(k+ 1)/’(t, y/k + 1) 

(96) 

(97) 
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r,p(T.£. (•/*+ 1 ) = 0 . fS3D. 



Proof: From (93* it follows that 
/>(T..r . ) 'k + I) ~ p(i*x.y/k) 

+£(t,a+i.x.a+ i)£[«-(A >\0 

f(k~ i)]B'(T,*+r'.u*+ 1) 

-B(r, A + l.x.A + \)E[r(k + \) 

*u'(t, >•/*)] - £[S(f,*A'Mk + I)] 

.£'(r, A +!,>•,*• )). 

£[5( T ,j:/*) f '(t+l)]=£[«(T.*A)J-(*/t)] 

•e:w’(A + 1) 

= L„(T,x/k)£:H'(k+l) 
and 

£[.r(A -«• l)fi'( T . V/A. )] = H(k + 1)lM;,(t. >•/*). 
Thus, we have 

p{ r.x, jr/A + I) = />(t. *, j'/A) + £(t, A + 1.x, A + 1 ) 
*r(A -+■ 1/A)£'(r. A + 1..V.A+ 1) 

— B{t, A + l.x.AH- l)tf(A + 1) 
•£.A.;,(t, y/A ) - I„(t, x/A ) 

■£'.H'(k + 1 )£'(t, A + l.y.A + 1). 
Substituting (79) into the above equation yields 
/>(t. x. y/k •*•' 1) = p(t. x. v/A) — x/A)i:jffi 

-Kk+ l)r “'{A + l/k)H(k +T)£.I;,(t. .v/A). 
In order to derive (96). note that from (81), 

. LjT'X>k)'l:=J(r.x/k-r 1/A*) 

and from the matrix inversion lemma (64), 

HIHPH- - Rf'H = (/ + H'R-'HP)~'H‘R-'H. 
Then we have 

H\k-r lir-'lA-f \/k)H(k M 1) 

— V ( A — l/k)R[k * 1) 
and 

p(~. x, y/k t- ))= p(r,x,y/k) - J(t. x/k + 1) 

•<t~\k + 1/A )R(k + l)y'(r, y/k + 1). • 

Multiplying each side of (94) by u(r, y/k + 1) and taking 
the expectation yields \\p{r, £. y/k 4 1) = 0, l £ dD. 
Thus, the proof of the theorem is complete, Q.E.D. 

Corolian A: J{ t, x/k ) satisfies the following relations. 

J(r, x/k+ ]) = A(7,x)J m {r+ 1/A; 4- 1) (98) 

and 


J(i + t. x/k) = Z)(r. x)J„(r/k) (99) 

where 


•U T /A) = 


J(r,x'/k) 

J(T.x m /k) 


A T. x) = p m ( T, X/t)£;^(t + J/ T ) 


(100) 

( 101 ) 


/>(r.x) = p„(T+ l,x/T)(p mm (T/T)£i) (102) 

Proof: Letting 4>(A 4 1) be given by <3>(A 4 1) = 
£i(/ + R{k + l)p mrn <* 4 from (88) and (89) it 

follows that 7 ( t, jc/A: 4 1) = p m { r, jc/t)$(t 4- 1 ) 4 >(t 4 
2) •.•*<* + !> and J w (r4 1/A4 1) = /v„(t4 !/t + 
1 ) 4 >(t 4- 2) • • ♦ 4>(A 4- 1). From the above equations and 
(75) we have 

J[r.x/k + 1 ) = p m (T,X/T)«I>(T+ 1 ) 


'Pmmi* + 1 A + 1)4.(T + V* + 0 


= */*)£•+( T+ V T )^ _, ( T+ V T ) 


+ 1A)-U T * 1 A + 0 

-A( 7 'X)J m { 7 + 1 /A- 4 - !)/s^ 

From (88) and (89) it follows that J{r 4 1. or A*) = 
p„,(r 4- 1. */r 4- 1 )&( t 4 2(* •* 0(A ) and JJr'k) = 
p mflJ (r/T)0(T 4 l)Otr 4 2) •** O(A). Thus, we have from 
the above equations 
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J ( t + 1 ,X/k) = p m (r+ J. Jt/r + 1)(/ * R(t + 1) 

•PmJ T + Ul))(p Mm (r/T)£'.)~'jjT/k) 
-p„{i + 

where the following equality derived from (66) has been 
used, 

P„( T+ 1. */ T+ 0 =/’»l( T+ >.*/' T ) 

•(/ + «(t+ 0^ b (t+1/t))"', (103) 

Thus, the proof of the corollary is complete. Q.E.D. 

Theorem )I: The optimal smoothing estimator is given 
by 


u(r % x/k) « 6 (t, */t) + 2 y(r.A'//)i!(/) 


l* T-* 1 

(104) 

r t C(r.f/A) = S(T,f). |6 0Z> 

(105) 

where 



(106) 


Furthermore, the optimal smoothing error covariance ma- 
trix function /?{x. a\ y 'k ) is given by 

p{ T. at. ,V/A) = />(t. X, .V,/t) - ^ 2 y(T. 3 

• e U/l~ \)AU)J\r>y/0 (107) 

r,p(T.(.>/A) = 0. (£30. (108) 

/*«w/ From (77) and (80). (104) can be directly ob- 

tained and from (96). (107) is clear. Thus, the proof of the 
theorem is complete. Q.E.D. 

VI, Summary of the Optimal Smoothing 
Estimators 

A Fixed ‘Point Smoother (x = fixed, k — r + 1, t + 

2. w * ) 

Theorem 12 The optimal fixed-point smoothing estima- 
tor is given b> 

u( t, x/k + ) ) = w( t, x/k ) + 7 ( T, x/k -f 1 )f (A + 1 ) 


(109) 

j(r.x/k+ i) = y(T,x/A')ir:V'(A*+ i/a) (no) 

+(A+1 k) = (l + R(k+\) Pm Jk+\))~' (111) 
J(r.x/i) =pjr,x/t) (112) 

r,i)(T.(/A + I) = s(r,f), (E3Z) (H3) 

r,J((.*/A+ 1 ) = 0 , (€3Z>. (114) 


Furthermore, the optimal fixed-point smoothing error co- 
variance matrix function p( t, x , y/A *f 1) is given by 

^(t. At, y/k + !)■•= p(t, x, y/k) -J(t, x/k + Ij^ Tj) 
.<(r; l/A)^(A + l)y'(r,>7A + 1) (115) 
r < p(T,(. 3 '/A + l) = 0, f£0i>. (116) 

5, Fixed-Interval Smoothing Estimator (A = /««/, t = A 

- 1, A - 2, • • • ) 

From Theorem 1 1 it follows that 
0 (t + 1 , X/k ) = fi( T + 1 .AT/T+ 1 ) 

+ 2 J(t 4 - l,x/l)i ! (l) ( 117 ) 

/=T + C 

and 

/>( t + 1 , j. y/A ) — p(x + 1* Jr, y/r + 1 ) 

/- t-: 

-J?(/)J'(t+ Ky/0- ( 118 ) 

Then the following theorem holds. 


IEEE: SYST n MAN, CYHERN 
Issue ‘ 

A 


■f*r- 


aT 



ORIGINAL PAGE ES 
OF POOR QUALITY 


Theorem 13 The optimal fixed-interval smoothing esti- 
mator is given by 

«(t + I ,x/k) — i5(t + 1, x/t + 1) 

•M(t+ 1. x)[u„,(t + 2/A' ) — 0,„(t + 2/t + l)] (119) 

r r fl(r+1,f//lt) = S(T +!.{)• %P‘ (120) 

Furthermore, the optimal fixed-interval smoothing error 
covanance matrix functior. ts given by 

P(t+ 1. x, y/k ) =/i(t+ 1. x, y/ t+ 1) 

-^(t+ \.x)(p mm (T+ \/k) 
~P*Ji + \/t))A'(t+ l,.r) 

(121) 

r,p(r+ hi. y/k) *0. (Ed D. (122) 

Proof. From (98) and (117) we have 
w( r + 1 , x/k ) = i5{ t + 1 , x/t + 1 ) 

k 

+4r^M 2 4,(t+1 //)?(/)• 

/e= T 4-2 

Bui from Theorem 1 1. 

w(t t* 2. x/k) = u(t + 2, x/t *f 2) 

k 

and from (43) and (59), 

w(t + 2. x/t + 2) = u(T*r2 t x/t + 1) 

+ F{r -r 2, x. t + 2)r( r + 2) 

= u(r ~ 2. x r + 1) 
hh*/(t — 2* x/t + 2)f (r *4* 2), 

Thus, wc have 

Utt2 tt 1) 

= 2 y m (Tf2 //)*</). 

/- T- 2 

Then we have 

w(t *+■ 1, x/fc) = «(t t 1,x/t+ 1) + A(t + 1.x) 

•Ifl«(r + 2/Ar) - tS m (r -+- 2/r+ 1)]. 
From (98) and (118). 

p ( t + 1 , x. y/k ) ~ p(r ‘ f 1.x. >*/t 1 ) 

—A(t + 1, x) i J m (r + 2/n0) 

^ (Ei*Z . — 7 

•W- 1)^(/R(t + 2//) 

•/<'(t+ !,>•)• 

From Theorem 11. 

^(t + 2. x. t r/A ) — /?(t + 2, x, j'/r + 2) 

- 2 i(f + 2,x//)rW/-WWT + 2 t )‘//) 

/=T*r 3 


6 / 
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and from (66). 

p{r + 2, A* v/t + 2) = /»(t + 2.x, >/t* 1) 

~ p m { r + 2. x/t l)^(f *r 2 /t <+ 1) 
•J?{T + 2)pi(T+2,>-/T+ 1). 

Taking into consideration that, from(103). /(t + 2, x/t + 
2) ss /> w (t + 2, x/t + 1)^(7 + 2/t + 1) and 

+ 2/If) - /> ww (t + 2/t *f 1) = 

- 2 J m (T + 2/m~\l/l-l)Ml)J:(r+2/l)' 

/°T-r 2 

we have />( t + 1, x, j'/A) = />(t + 1, x. y / t + 1) - ^(t 
+ J. xX A h M (t + 2/A ) - />*,*,{ t 4* 2/t + 1)]/T(t + h .v). 
Since the boundary conditions (120) and (122) are clear 
from (105) and (108), respectively the proof of the theorem 
is complete. Q.ED. 

C Fixed-lag Smoothing Estimator (t = k + 1. A = A + 1 
4* A, A » //xed ) 

From Theorem 1 1 we have 
u{k + 1, x/k + 1 *+■ A) - ri(A + 1, x/A + 1) 

i-rl+A 

*.««-■ + 2 J(A 4- 1 , X//)r(/) (123) 

/« A — 2 

*p(A 4* 1, x, v /A **■ 1 + A) 

= p(k + lx. v/k + 1) 
a-i-a 

- 2 + 1 

•*(/./- l )«(/)/( A + i. jr//). 

( 124 ) 
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Then the following theorem holds. 

Theorem 14. The optimal fixed-lag smoothing estimator 
is given by 

u(k + \,x/k + 1 + A) 

= £,u(k.x/k + A) + C(x,k + 1,4) 

■F m (k + I + A/* + 1 + AM* + I + A) 

+ QJk,x)(p m „(k>W,)~' 

■(u„(k/k + A)-u m (k/k)) (125) 

r t 4(A. \.t'k+ i + A) = s(A + i.{). fea d 

( 126 ) 

where 

C(x. k + \, A) = A(k + Ux)Ajk + 1). 

••• ,A m (k + A ) (127) 


A m (k) = 


Afk.x') 

[^(/.jc m ) 


Furthermore, the optimal futed-lag smoothing error covari- 
ance main); function p{k 4 1 t x,y/k+ 1 4 A) is given 
bv 

p{k 4 1, x, \/k 4 1 4 A) 

= p{k+ \.x.y/k)~ C{x,k+ 1. A) 

• F„( k 4 ] -T A fk 4 1 4 A )H(k 4 1 4 A) 

•p K Jk - 1 A /k 4 A)C'(,v» k 4 1. A) - D{k, x) 


■[p m Jk'k)-p m Jk'k + A))D-(k.y) (128) 

T t p(kt \.i.y/k + 1 + A)= 0. f £ 32). 

(129) 


Proof From (43) and (59) we have 
u[k - 1..V/A + 1) = e,uU,jc/A) 

+y(A + 1 . jc/A - \)v(k-r 1). 
From (123) and the above equation it follows that 
u(k + 1 ,x/k+ 1 + A) = £.«(*, x/k) 

t-a 

+ 2 Ak + \,x/k)f(l} 

/= A 1 

»(A-+ l,x//+ 1 + A)P(A + 1 + A). 
From (88) it follows that 
J(k + 1. x/k + 1 + A) = />„,(* + 1. x/A + 1) 

■i',t(k + 2/k+ 1) 


•£;*(* + 3/A + 2) ••• 


•e;^(A- + 2 + A/A-M +A). 



4 / 


IEEE: SYST., MAN, CYBEBN. 
Vol. ,SSU8 

Au: CHUii. l 6il 

Galley No, ( 7 


ORIGINAL PAGE 13 

Substituting (75) into the right side yields OF POOR QUALITY 

J(k + 1.x/* + 1 4 A) 
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/(A: 4 \,x/k 4 1 4 A) 

>•(* 4 1 + A) = C(x. At 4 l. A )p mm (k + 1 + A/A 4 1 4 A). 

(130) 


From (99) it follows that 
2? 3(*+l.x//)r|/)= *2 Fm(* + l.x/At) 

/*= A — t 

(^JA/A)£;)"'7 m (A-//)0(/). •>/ 

But from (29) we have 

/>„(* 1 . x/k ) = i,pjk,x/k)i: 4 QJk.x). 

From (9S) and (99) we have 

Ak.x/1) 

-A[k.x)Jjk - 1//) -p„(r, x/r)£’,( p mm {k/k)i'.}~ t J m (k/I). 

Then it follows that 

2 Ak - i.x//)f(/) = e, 2 Ak.x/nnn 

M-l (=A-I 

■rQJk.x) k 2 (p m Jk k)i‘.)~'j m (k/l)HI ) 



and 
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A(A + \,x/k + J + A) 

~£,Ci{k,x/k+ A) + C(x,A + l.A) 

+ I + A/A + 1 + A)i*(A + I + A) 

+$Jk'X)(p m Jk/k)£'.)~' *2 JJLk/nm- 
(**1+ \ 

Bui from Theorem 1 1 we have 

k + X 

1 

Thus we have (125). From (65) and (124) it follows that 
p[k+\.x.y^*pk)«p(k*Ux.y/k)-J,^J 3 i 4 / 
where 

/,= V M+UX/W~'V/I- l)*U) 

•J'(k + 1..V//) 

J : = J(k + \.x/k+ 1 +AH“'(*+ 1 + A./A + A) 

•/?(A + ) + A)/'(A + 1. ,v/A + 1 + A). 

From (75) and (130) we have 

J ; = C(x, k + 1. A )p„Jk + 1 + A/A + 1 + A) 

■R(k 1 + A)p mm ( A -r 1 + A/A + A) 

•C'(.t,A - l.A) 

= C(x, A + 1. A)F m (A + 1 + A/A + I + A) 

■//(A + 1 ■+ &)p mm (k 1 + A/A + A) 

>C\y,k+ l.A). 

Substituting (99) into 7, yields 

/, = D(k.x) k j? JJk/l)<T'(ri- 1) 

/=*- l 

*K(/)J;(A 7}2T(*. >). 

But from Theorem 1 1 we have 

p{ A , x. v/A + A ) - p( ft . *, r/A ) 

= - *£* y(A,jc//K -'(///- i)*(/).n*oy7) 

I 

and 

p„„{k/k * A) - p mm (k'k) 

= - *2 Jjk/IH-'V/I- \)R(IK(k/l). 

/= A- 1 

Then we have 

/>(A ■*+ 1, y/A + 14 5) 

“ p( A +• 1* x, .v/A ) - C(x, A *+ 1, A) 

•F m (k + 1 + A/A + 1 + A)//(A + 1 + A) 

'Pmm( A + 1 + A/A + A )C‘(y, A + 1. A) - Z>(A,x) 
•[p„„(k/k)-p mm {k/k + L)]D’{ k,y). 

Since the boundary conditions (126) and (129) are clear 
from (105) and (108). respectively, the proof of the theorem 
is complete. Q.E.D. 

Kelly and Anderson [18] proved that the fixed-lag 
smoothing algorithm of Theorem 14 may be unstable, but 
Chirarat tan anon and Anderson [19] derived a stable ver- 
sion of the algorithm. It is possible *o derive a comparable 
version here, although stability problems should not arise 
in our use of the algorithm of Theorem 14 as long as it is 
used over a finite time interval. 
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IX. Application to Estimation of Air 
Pollution 

Distributed parameter estimation theory has recently 
been applied to simulated air pollution data to demon* 
strate the capability of estimating atmospheric concentra* 
tion levels from routine monitoring data [10], 111). A 
problem identified in thesej early studies was how to 
specify the statistical properties of the assumed system and 
observation noise. In this section we expand upon the prior 
studies in two respects. First, we consider actual monitor* 
ing data for sulfur dioxide (SO,)* in particular those mea- 
sured each hour during the period December 1-31, 1975 at 
four locations in Tokoshima Prefecture. Japan (sec Fig. 1). 
Second, wc applv the method of Sage and Husa [12] to 
estimate the unknown noise covariances in the system 
equation and measurements. 

Hourly sulfur dioxide data are available at the four 
locations shown in Fig 1 for the period December 1-31, 
1975, The data for day k at location i may be denoted by 
e k {x\t). It is useful to average the data for December 
1-30 to produce 

w*'.o>4 1 f, (**.») 03i) 

M A*1 

where we will consider December 31 as a day to test the 
algorithms. 

If it can be assumed that the wind flows are such that 
there are no north-south variations of conceniration and 
that vertical mixing is rapid enough to eliminate variations 
of concentration with altitude, then the region can be 
considered to be one-dimensional along the east-west coor- 
dinate. The S0 2 concentration at any particular time can 
be assumed to be described by the atmospheric diffusion 
equation [13]. 

where f is the wind velocit), a is a diffusion coefficient, 
and 5 is the rate of emission of SO* as a function of 
location and time. 

Equation holds at any instant of lime, but we 
desire an equation governing the monthly mean concen- 
tration <c). Although no such equation exists, we can 
formally average (132) over the 30 realizations (days) to 
produce 

* - (133) 


3 <£> 

3/ 


■(*>-(• 


3* 3 


+ S. 




3x 


= 0. x = 0, 1 


f-O. 

3x 
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Fig. 1 


One object will be to estimate the diffusion parameter a. 
This parameter will in general vary with location and time 
of da>* although for simplicity we seek a constant value for 
the month. Thus, the first term on the right side of (133) 
becomes a 3 2 (u)/3.r 2 . We can form the residuals, u — c — 
<c> and z = e — (e), By subtracting (133) from (132) we 
obtain 

3u, .3c /..ScV 9 2 u /, « A \ 

37 + ^-( f 3?) = fl i7- (134) 

Since wind data are not available with which to evaluate 
the second and third terms on the left side of (134) let us 
rewrite (134) as 

3u 


(135) 


where w(x, /) includes those unknown features associated 
with the velocity terms. 

The boundary conditions on (132) are 

3c 


(136) 


expressing the assumption that there is no diffusive flux of 
SO* into or out of the region at the boundaries. After 
averaging and forming the residual, (136) becomes 


*/ 


x = 0. 1. 


(137) 



A fc fe 


(* tl llM 


The problem is now to estimate u(x> r) based on the 
data, 


:(x'.t) = u{x\t) + c,(t). ( « 1,2, 3*4. 0 38) 


^riirtlWikL PAGE 521 


Since hourly data arc available. (135) can be cast ime* the 
discrete-time form (1), 

u(k + 1.x) = £ a u(A\ x) + H’(jc* /) (139) 

with i M =s I + fl3 : /3x 2 . Observation error is estimated 
from the mean square error of predicted values and ob- 
served data. 

M{t)n± 2 (.-,(*) -<,( k/k - !)) 5 . f* 1, 2.3,4 

** k *» I 

(140) 

An index ol overall estimation error is 

J* 2 w<0. 04i) 

To apply discrete-time distributed parameter estimation 
theory to predict air pollution levels, we must consider 
three problems. The first problem is how to simulate the 
distributed parameter system. The second is how to de- 
termine the covariances of system and observation noise. 
The last is how to determine the diffusion coefficient a. For 
the first problem we use the Founcr expansion method and 
approximate the original distributed parameter system by a 
finite-dimensional system. For the second problem, we 
apply the algorithm of Sage and Husa [12] that necessitates 
the simultaneous application of the optimal filtering and 
smoothing algorithms For the third problem we apply the 
maximum likelihood approach in the smoothing form [14). 
We now consider these problems in more detail. 

Fourier Expansion Method It is well-known that the 
state u(k*x) of the distributed parameter system (139) 
with boundar> condition (137) can be represented by using 
the eigenfunctions £,{a ) as follows, 
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u(k.x) = 2 ",(*)$.(*) 

4= 1 

(142) 

where ! 


£,$,(*) = x £i(0. 1) 

1 1 

T'»' <-M 


and 

(143) 

Au/Ed: placement of eq. no? 

J\lx)<t>j[x)dx = S tJ . 



A, is the eigenvalue of £ x corresponding to </>,(*). In this 
case, it is easily seen that the eigenfunction <f>,(x) and the 
eigenvalue are given by 

<fp,(.x) = l, $,(*) = »£ cosr.K, i = 2y. 


and 

X, = 1 — or s (» — l) s , /= 1,2.--*. ("**) 

Then «(t, x/A*)* P(t. xf /A), and A(r, x) can be rep- 
resented as follows: * 



u(;. x/k) = 2 &,j(r/k)$,(x) 

I* I 

p(T.x,y/k)= 2 P,j(r/k)<t>,(x)4>J(y) 

.,y-i 

A(t.x) = 2 a t (T)4.,(A-)- (* 45 ) ^ / 



Let us approximate these infinite expansions b> the first A' 
terms and define the following matrices and vectors. 

fitT/A) =* Col[a,(r/A*).***-fiA(r/Jt)]. 

A(r) a Co|[<3,(t).-**.o v (t)]. 

A = diag[A,.-->, X A J. 



‘ Pn(r/k).' 

••‘‘P i.v(t 

F(t/A) = 

J 

* 


_Ps\(' T / k ')-' 

•••'Pssi' 


4h (*)."* 

. 9tx( A) 

Q(k) = 




4s ,(*).••• 

.?a.v(A). 


’ *,(*').••• 

.♦at*') ' 

4> = 

: 






and 


where q,,(k) denotes the (<, ylth Fourier coefficient of 

Q{k.x,yY 

Then, from Theorems 3-5 we have 
u(k + \/k + 1) = Au(A/A) + F(k + l)r(A + 1) 

F(A + 1) = P(k + 1/A)4>7/'(A + I) 

[tf(A + 1 )4>F(A 


A + \/ kW@) 


■{k ~ 1) + R(k + 1)]* 

P{k + I/A ) = AP[k / k).\' + Q(k), 

P(k + 1/A + 1) = (/ - F(k + l)f/(A + 1)<J>) 

■P(k* 1/A). (146) 

Furthermore, from Theorem ]2 we have 
u(r + 1/A) — $(t f 1/4-1) + A( r + 1)4 >(u(t + 2 /A) 
— u(t 2 t + 1)). 


ORIGINAL PAGE 1$ 
OF POOR QUALITY 



IEEE: SYST., MAH, CYBERN. 
VoU Issue 

* 

Aoi: C)/y*CL t\s 

Galley No. tD Z>- 


13 

x / 


A[ r+ 1) = F(t — 1 /t — 1)A?-'(t+ 1/t)4>-', 

F(t+ 1/A) = P(t -*• l/r+ 1) 

-A{r + 1)4>(F(t-t 1/A) 

-?{t+ i/ t ))4.^'(t+ 1), (i«n) 

Note that the fixed-interval smoothing estimator does not 
depend on the matrix which reflects the effect of sensor 
location. 

Determination of the Noise Covariances: In order to 
determine the unknown covariance matrices of the system 
and observation noises, we adopt Sage and Husa’s algo- 
rithm (12) given by 

d(*)»T 2 (<JUA)-Ai0 -1/A)) 

’(*(j/k)-Au(j-\/k)Y (148) 

and 

*(*)=T 2 (:U)“HU)*0U/*)) 

“ J * I 

•(*(;> 049) 

where @(k) and R(k) denoted the estimated values of 

Q(k) and R(k), respectively. Note that in these identifier "tku, 

lion algorithm the fixed-interval smoothing estimate u(j/k ) 

is used. 

Identification of the unknown Parameter a: To determine \J 

the unknown parameter a we use the maximum likelihood 
approach in smoothing form (14). The log-likelihood func- 
tion y(Ar; a) is given from [14] by 

y(k\ a\— 2 “ (Tbui 


(150) 



where 

w= -If In (2* I- i In del ~ I .a) 
r<**~ - 2 {p'UitfJ/rVt;;*) * (u(j 'k >a) 

j a * 

-tfO- 


•<*0 *>*)-«;- 1/A.fl))} 

iij.c)z* :{j)~ Hij^ij j~ Uoi $ 

£,(./ 7 - 1. 0) « a)i*'0: <?)]. 

where p is the dimension of ;(A). and u(j k ~ her) 
denote* u< ; V. - 1 j under the condition that the unknown 
parameter is assumed to be a 
To maximize ytA.c) we use the following gradient 
method 


0 * I *= Q, •** C?U)Y,(A; tJ, ) 
0y(A.0) 


y f {A,0,) - 




051 ) 


where Cm is a suitable matrix Therefore, we adopt the 
following recursive algorithm to idennfs the unknown 
parameters Q, R, and a 

1) Make an initial guess a ( , of a, 

2) Compute Q »u ( 1 and R\ a K t bv using (14S) and (149) 

3) Compute c b> using [\ 50) 

4) compute ( [ha 1 and A{a t ) bv using (14S) and (149). 

5) Reiurnitbree b> changing mo / - 1 and repeat until 
these values do not change 


Suberic a! Retain We use the observed data from De- 
cember 1-30 to identiH the unknown parameter a and 
noise covariances Q and R After four aerations the algo- 
rithm for determining a converged to the value a = 0.001 
The Fourier expansion hat been truncated at V s 4 The 
estimated diagonal elements of noise covariance matrices 
are 
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Q t: - 6,44 A,, = 0.29 

Q :2 =}A0 R :: ~ 0.61 

qZ = $.75 rZ= 1.96 

Q m = 2M Rt 4=1-34. 

To consider the effect of the nuipjer and location of 
monitoring stations, we assume that we have data at oniv 
one monitoring station. In this case from the previous 
results of Kumar and Semfeld (15) and Omatu ei ah (16) 
we expect that the optimal sensor location is closest to the 
boundary Thus, cither x 1 or x 4 is the optimal single sensor 
location among the four monitoring stations, xk x 2 < jc\ 
In Table 1 we show the values of PA(i) and J for several 
moniionng stations We see that Aizurru or Matsushige is 
optimal for the one-point sensor location case. Similar 
conclusions hold for two or three monitoring stations. 
Fin allv. wc illustrate the actual observation data and one- 
hour ahead predicted values for December 31 in Figs, 2-5 
for Aizumi. Kjtaji^ma, Kawauchi. and Matsushige. respec- 
tively. 

Comparison wth Other Approaches; It is of interest to 
compare results of the present filtering and smoothing 
approaches with others available for air pollution estima- 
tion. We consider, therefore, the same S0 2 estimation 
problem by the following methods. 1) XR-modek 2) per- 
sistence, and 3) weighted ensemble, 

The^R-model method is based on the following AR{p) 
model 


A/ 


Tabic I, Figs. 2-5 


u l" = + a i u l'-i + + V*-/> + 

/= 1. 2.3.4 (152) 


n 



where the are the concentration levels at time k and 
at monitoring station x'. a v a t are the corre- 
sponding /^parameters, and the e; "& arc re«tduab V>e 
used the Levinson algorithm to determine the AR* 
parameters while the optima! ord ** p of the ^/Nprocess is 
determined by using the minimum Akaike’5 information 
entenon {AlCi |20J Then the one-hour ahead predicted 
concentration is given b> 


* ojul'ii * 

and the prediction error variance is 


J ~ 


4 ( 


74 

2 (C- 


(153) 

(154) 


. A M « M I 


t/y 


ICEC: $YST,» MAN, CY3£flN. 
Vol, ‘ issue / 

Au; ^ 1 

Galley No 


Tab! II shows* the .4 ^parameters and minimum AJC 
value at ca*h monitoring station Table 11 

The persistence method consists mercl> of using the 
observation data as the one-hour ahead prediction 
value u” // 

The weighted ensemble method uses the mean of the 
past observation data at each time k weighted by a linear 
function of the source strength as the prediction value at 
lime k Ba^ed or. the number of emission sources, the 
weighting functions are assumed here to be 0 15* 04). 0.26, 
and 0 IS at x\ x x\x\ respectively Table III shows the 
performance criteria of the four methods From Table III 
we can see that the present method possesse> almost the 
same accuracy as the AR-model method B> multiplying 
each eigenfunction coefficient by the corresponding eigen* 
function 3nd summing them, however, the present method 
enable 1 ' us to estimate concentration* over ;&« enure re* 
gr*n Therefore, the present method is more poweWu! than 
the A /?- model method Table III 


IX Cqvcu $io\s 

Optimal Citimaiors for discreie-une distributed parame- 
ter ss stems have beer, derived based on Wiener- Hopf 
theory A notable point of the present work is that the 
smoothing estimators have been derived by the same ap- 
proach as the filter, thus providing a unified approach for 
thi* clas^ of distributed parameter estimation problems. 
The estimation algorithms have been applied to the prob- 
lem of predicting atmospheric sulfur dioxide levels in the 
Tokushima prefecture of Japan 
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Fig 1 Map of Tokushima Prefecture Japan The four air pollution 
monitoring nation showm arc located as follows a* Amirm. 
Kitaima Kattaueh; x* Matsushige Sources of sulfur dioxide 
has c been lumped according to the three sources sues indicated by the 
open circles p-v 30*50 mV h 10-30 nvvb • < 10 m'/b 

v^> v 

Fig 2 Measured and estimated sulfur dioxide concentrations on De. 
ccmber 31. 19?5 at Auuroi monitoring station (a ’) 

Fig 3 Measured and estimated sulfur dioxide concentrauons on De- 
cember 31. 19'i at Kitajima monitoring station l a 3 ) 

Fig 4 Measured and estimated sulfur dioxide concentrations on De* 
ccmlet 31, 1975 atKittauchi momtonng station (a 3 I 

Fig 5 Measured and estimated sulfur dioxide concentrations on Dc» 
cember 31. 1975 at Matsushige monitoring station (aV| 
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ESTIMATION OF ATMOSPHERIC SPECIES CONCENTRATIONS 
FROM REMOTE SENSING DATA + 
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ABSTRACT 

A basic problem in the interpretation of atmospheric remote sensing data 
is to estimate species concentration distributions. Typical remote sensing 
data involve a field of view that moves across the region and represent inte- 
grated species burdens from the ground to the altitude of the instrument. 

The estimation problem arising from this special measurement configuration is 
solved bas.id on the partial differential equation for atmospheric diffusion 
and Wiener-Hopf theory. The estimation of the concentration distribution 
downwind of a hypothetical continuous, ground-level source of pollutants is 
studied numerically. 
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I, Introduction 

In the remote sensing of atmospheric species, a ground-, aircraft-, or 
satellite-based platform scans a region of the atmosphere and measures the 
species burden within the field of view. An object of atmospheric remote sen- 
sing is to reconstruct species concentration distributions over a region based 
on the data available from the instrument. 

There exist two recent studies that assess the capabilities of remote 
sensing for monitoring regional air pollution episdoes [1,2]. Diamonte et al . 

[3] developed theoretical results for the estimation of point source plume dis- 
persion parameters from remote sensing data. In a similar vein, Kibbler and 
Suttles [41 studied the estimation of unknown parameters in a pollutant disper- 
sion model by comparing model predictions with remotely sensed data. No results 
have yet been reported in which actual remote sensing data have been used to 
estimate species concentration distributions. 

The present paper deals with the theoretical foundation of estimating atmo- 
spheric concentration distributions from remote sensing data. Since the atmosphere 
is a three-dimensional system, mathematical models of pollutant behavior are of the 
distributed parameter type [5], Remote sensing data usually represent spatial 
averages of concentrations, so that the estimation problem concerns a distribu- 
ted parameter system with spatially integrated, scanning data. Although dis- 
tributed parameter state estimation has been considered extensively (see, for 
example, [6] and [7]), such problems with scanning and spatially integrated 
measurements have not been considered previously. The purpose of the present 
paper is to derive the required optimal estimators for the scanning and spatially 
integrated measurement case by a unified method based on the Wiener-Hopf theory. 
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In Section II, we define the remote sensing data analysis problem mathema- 
tically. Sections III-VI are devoted to derivation of the optimal prediction, 
filtering and smoothing algorithms for the problem by Weiner-Hopf theory. Fi- 
nally, in Section VII we present a detailed numerical example of estimating the 
concentration distribution downwind of a continuous, ground-level line source to 
illustrate the application of the thoery. 

II. Problem Statement 

We consider a single atmospheric species (nonreactive), the mean concen- 
tration u(t,Xj»X 2 ,x 3 ) of which over a certain region is described by the follow- 

"k 

ing form of the atmospheric diffusion equation [5], 


3u ^ w 3u ^ y 9u _ 9 
9t V 1 9x 1 v 2 3Xg ~ 9X.J 




+ w(t,Xj,x 2 ,x 3 ) 


( 1 ) 


where and V 2 are the mean velocities in the x^- and x^-directions , respec- 
tively, K y ( x^ ) is the vertical turbulent eddy diffusivity, and w(t,x 1 ,x 2 ,x 3 ) is 
a random disturbance accounting for inaccuracies inherent in the basic model. 
The initial condition for (1) is u(t o ,x^,X 2 >x 3 ) = u q (Xj,X 2 »x 3 ) , and typical 
boundary conditions are 


" k v (x 3^ ax 3 = ^(t,x 1 ,x 2 )> 


9u 

9x, 


= 0, 


x 3 = 0 
x 3 = h 


( 2 ) 


where S^Xj.Xg) is the ground-level species source emission rate, presumably 
a known function, and h denotes the upper vertical boundary of the pollutant- 
containing region, for example, the base of an inversion (stable) layer. For 
convenience, we denote the coordinate vector by x and let 


L [•] = 
x l J 


v Mil 

1 3x 1 


2 3x 


iLJL + 


9x. 





k 

In this form of the atmospheric diffusion equation, turbulent diffusion in the 
horizontal direction is neglected relative to transport by the mean flow, a 
common assumption in treating atmospheric diffusion problems [53. 
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Assume that the concentration of a species over a fixed spatial domain D 
with its boundary 3D is of interest. Let us define the operator r^, 5 £ 3D 
as follows, 




■ w ft 

SI-1 
3X 0 * 


a[*i 


x 3 ‘ 0 
*3 ’ h ’ 


Let S(t,E) be 


S(t,£) 



Thus, (1) can be represented as 


= 0 
= h. 


— ■ 1" -^ = L u(t,x) + w(t,x) (3) 

0 A A 

and (2) can be written as 

r^u(t,C) = S(t,?), ^ £ 3D. (4) 

We assume that the initial condition u (x) can be represented as a Gaussian 
process with statistics, 

E[u 0 (x)] = u Q (x) 

E[(u 0 (x) - u Q (x)) (u Q (y) - u Q (y)j ] = P Q (x,y) (5) 

and the random disturbance w(t,x) is stochastically independent of u Q (x) and is 
a white Gaussian process with statistics, 
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E[w(t,x)] =0 

( 6 ) 

E[w(t,x)w(s,y)] = Q(t,x,y)5(t-s). 

We assume that the remote sensing measurements are taken at time t^ over 
a view volume D(k) consisting of M pixels, as shown in Fig. 1. Since the sens- 
ing platform may be in motion, the field of view, in general, moves with time 
across the entire spatial domain D. We assume that the shape and extent of the 
field of view D(k) remain fixed and only the location of the centroid of each 
pixel changes with time. The ground-level location of the centroid of each 
pixel of D(k) is denoted as (x^ k \ x^^.O), m = 1,2,..., M 

We are interested in considering the vertically integrated measurement 
given by 


Z m(k) (t k’ n ) = / n ^m(k)^ x 3^ u(t k ,x T (k) ’ x 2 (k),x 3 ) dx - 


v(t k ,xf k) ,xj (k) ,h n ) 


(7) 


m = 1,2,..., M, n = 1,2,..., N, hj < h 2 < ... < h^ 
k = 1,2,... 

where 1s an altitude-dependent instrument weighting function, and h n is the 

vertical position of the scanning sensor. Physically, ^^(t^n) represents 
the vertically-integrated species concentrations within each of the M pixels, 
indicated by m(k), at each time, t^, from an altitude of h n - v(t k >x^ k ^ ,x^ k ^ , 
h n ) represents measurement errors. 


-5- 


Some comments concerning the measurement configuration shown in Fig. 1 
are in order. Ordinarily remote sensing from an airborne platform would be 
carried out at a single altitude. In such a case, it is not possible to esti- 
mate the concentration distribution between the platform and the ground based 
only on the integral of the concentration. Sakawa [8] and Koda and Seinfeld 
[9] have shown that in problems of this nature it is impossible to estimate 
the state uniquely based on integrated measurements from only a single sensor 
position since the required distributed parameter observability condition does 
not hold. Therefore, the estimation of species concentration distributions 
necessitates traverses over the region at different altitudes. From a practi- 
cal point of view this requirement restricts this type of monitoring to air- 
craft platforms, which, for purposes of measuring air pollution, are the most 
useful. Considering that atmospheric concentration distributions change gradu- 
ally and that airplane speeds are fast, the configuration sketched in Fig. 1 
implies that repeated measurements at several altitudes are possible using only 
one :irborne platform. 


~ 6 ~ 

In order to represent (7) more compactly we introduce the following 
notation: 


m(k) " 

d n (t k ,x 3 ) 


ij 



d(t k ,X 3 ) 


Z(t k ,n) = 


Z(t k ) = 


v(t k> n) = 


x(m(k)) = (xf k) , 
^m(k)^ x 3^ x 3 ~ h n 

w 0 » x 3 > ^n 

^l(k) (x 3^ * * ? 

• * . 

» 9 * 

* ‘ * J M(k)^ X 3^ 

u(t k , x(l (k) ) , x 3 ) 

u(t k , x(2(k)) , x 3 ) 

» 

f 

u(t k , x(M(k))»x 3 ) 

J (t k> x 3 ) 

^J N (t k ,X 3 )_ 

Z l(k) ( V n) 

Z M(k)^k’ n ^ 

r z (t k .D~ 


Z(t k> N)J 

v(t k , x(l (k) ) ,, h n ) 

» 

v(t k , x(M(k) ) , h n ) 
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v(t k> = 


v(t k ,l)~' 

* 

v(t k ,N) 


Then (7) can be represented compactly as 



( 8 ) 


We assume that v(t k ) is independent of w(t,x) and u Q (x) and is a white Gaussian 
process with statistics, E[v(t k )] = 0 and E [v(t k ) v"( t^) ] = R(t k )5 k ^, where " 
denotes the transpose operator and R(t k ) is an MNxMN positive-definite matrix. 

The problem considered here is to estimate u(t,x) over D on the basis of 
the measurement Z ( t^ ) , o = 0, 1,..., k. The novel aspect of this problem from 
the point of view of distributed parameter estimation arises because of the 
scanning and vertically integrated nature of the measurements. In what fol- 
lows, we use k instead of t, as long as there is no ambiguity. 
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III. Estimation Problems and Wiener-Hopf Theory 

Let us denote the estimate of u(t T> x) based on the observation data 
Z(t ), a = 0,1 k by u(t T ,x/t k ) which is given by the following linear 

transformation of Z(t ), a = 0,1,..., k, 

0 

k 

o(t T ,x/t k ) = y F(t t ,x,t 0 )Z(t a ) (9) 

o=0 

where F(t ,x 5 t ) is an unknown MN-dimensional row vector called the estimation 

T O' 

kernel function. When there is no ambiguity, we write (9) compactly as 

k 

u(t, x/k) = Ya F(t ,x,t )Z(t ). (10) 

a=0 too 

Furthermore, we denote the estimation error and error covariance functions by 
u(t T ,x/t k ) and P(t T ,x,y/t k ) , respectively, where u(t T ,x/t k ) = u(t T ,x) - u{i T ,x/t k 
and P(t t ,x,y/t k ) = E[u(t T x/t k )u(t T ,y/t k ) ] . The estimate u(t T ,x/t k ) that mini- 
mizes J(u) = E[u(t T ,x/t k ) 2 ] is said to be optimal. Note that by using 
p (t T ,x,y/t k ), J(u) can be rewritten as 0 ( u) = P(t T ,x,x/t k ) . 

To clarify the differences between the prediction, filtering, and smooth- 
ing problems, we express F(t ,x,t ) differently for each problem as follows: 

T L. 

(i) Prediction (t > t k ) 

k 

u(t,x/t R ) = Y* A(t,x,t a )Z(t 0 ). 


(ID 
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(11) Filtering (t T = t^} 

k 

u(t k ,x/t k ) = E F(t k ,x,t ff )Z(t c ). (12) 

(iii) Smoothing (t T < t^) 

k 

u(t T ,x/t k ) = 52 B(t T ,t k ,x,t 0 )Z(t 0 ). (13) 

Here we use three temporal arguments t , t k and t 0 for the smoothing kernel 
B(t^ ,t k ,x,t c ) since these parameters should be changed according to the measure 
ment data acquisition time. Then the following theorem can be proved similarly 
to that of [6] for the continuous-time observation case. 

[Theorem 1] (Wiener-Hopf Theorem) 

A necessary and sufficient condition for the estimate u(t T ,x/t k ) to be optimal 
is that the following Wiener-Hopf equation holds for £ = 0,1*..., k an-’ 
x t D = D . 8D, 

k 

52 F(t ,x,t )E[Z(t)Z"(t )) * E[u(t 4 x)Z"(t )] , (14) 

o=0 T 0 o A T 4 

or equivalently, for e = 0,1,..., k and x € 0, 

E[0(t T ,x/t k )Z"(t ? )] = 0. (15) 

[Corollary 1] (Orthogonal projection lemma) 

The orthogonality condition, E[Q(t T# x/t k )u(t ,y/t k ) ] = 0, x,y € B, holds where 

t is any time instant such as t < t, , t = t. . or t > t, . 
n J nk nk n k 
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[Proof] Multiplying each side of (15) by F^t ,y,t r ) and summing from 

n h 

x, - 0 to x, = k yields 

k 

E[0(t T ,x/t k ) £ Z'(t ? )F'(t ri ,y,t ? )J = 0. 

C— 0 

Using (9) in the above equation yields the desired relation completing the 
proof of the corollary. Q.E.D. 

[Lemma 1] (Uniqueness of the optimal kernel) 

Let F(t ,x,t ) be the optimal kernel function satisfying the Wiener-Hopf 
equation (14) and let F(t ,x,t ) + F A (t T ,x,t a ) be also the optimal kernel func- 
tion satisfying the Wiener-Hopf equation (14). Then it follows that F^(t t ,x,t CT ) 
c = 0,1,..., k and x£ D, i.e, the optimal kernel function is unique. 

In order to consider the prediction, filtering, and smoothing problems, 
separately, we rewrite (14) using the notation of (11) - (13). 

[Corollary 2] The Wiener-Hopf equation (14) is rewritten for the prediction, 
filtering, end smoothing problems as follows: 

(i ) Prediction (t > t^) 

k 

52 A(t,x,t )E[Z(t )Z'(tJ] = E[u(t,x)Z"(t )] (16) 

o=0 o a q x, 

for x, = 0,1 ,. . . , k and x £ D. 

(ii) Filtering (t T = t^) 

k 

52 F(t k ,x,t 0 )EtZ(t 0 )Z'(t ? )] = E[u(t k ,x)Z'(t c )] 


(17) 
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for t, = 0,1, , . . , k and xt D, 


(iii) Smoothing (t T < t k ) 
k 

YL B(t T ,t k ,x,t ff )E[Z(t 0 )2"(t 5 )1 = E[u(t T ,x)Z"(t 5 )] (18) 


for x, - 0,1 , ... , 


k and x € D, 
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IV. Derivation of the Optimal Prediction Estimator 

In this section we derive the optimal prediction estimator by using the 
Wiener-Hopf theory in the previous section. 

[Theorem 2] The optimal prediction estimator is given by 


3u(t,x/tJ 

_____ = L x u(t,x/t k ), t> \ 

r_u(wt k ) = s(t,c) * c e 3D. 


(19) 

( 20 ) 


[Proof] Differentiating (16) with respect to t and substituting (3) yields 
k 


E 

o=0 


3A(t,x,t ) 


at 


E[Z(tJZ'(t • = L £[u(t,x)Z"(t ) ] 

U A L, 


where the independence of w(t,x) and Z(t ) is used. Substituting (16) into 
the above equation yields 


Z F A (t,x,t a )E[Z(t a )Z'(t ; )] = 0 


where 


F A (t.x,t a ) = 


3A(t,x,t ) 

§t L x A(t,x,t a ) 


From Lemma 1 we have 


3A(t,x,t ) 

= L A(t,x,t ). 


at 


( 21 ) 
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Differentiating (11) with respect to t and substituting (21) yields (19). 
Since the forms of and S(t,£) are known, the predicted estimate u(t,x/t k ) 
also satisfies the same boundary condition (4). Q.E.D. 

[Theorem 3] The optimal prediction error covariance function P(t,x,y/t k ) is 
governed by 

3P(t,x,y/t.) 

9t — ~ = (L x + i - y ) p ( t > x »y/ t k ) + Q(t>x.y)» 

r ? P(t,?,y/t k ) =0, K € 3D. 

[Proof] From (3), and (19) we have 
3u(t,x/t.) 

^ — = L x Q(t,x/t k ) + w(t,x) (24) 

and from (4) , and (20) 

iyj(t,£/t k ) = 0, K € 3D. (25) 

Differentiating the definition of P with respect to t and using (24) yields 

3P(t,x,y/t k ) v 

^ = (L x + L y )P(t,x,y/t k ) + Z(t,x,y) 

where 

2(t,x,y) = E[w(t,x)u(t,y/t k )] + E[u(t,x/t k )w(t,y)] . 

Let the fundamental solution of L be G(t,a,x,y), where 

A 


( 22 ) 

(23) 
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B l x 6(t,c,x,y), 

r ? G(t,o,c»y) «* s(t,0, 3D 

G(c,r,x,y) = 6(x-y). 


Then u(t,x/t k ) of (24) can be represented in terms of G(t, o,x,y) as follows, 



Substituting (26) into Z(t,x,y) and using (6) yields E(t,x,y) = Q(t,x,y). 
Multiplying each side of (25) by u( t.y/t^) and taking the expectation yields 
(23). Q.E.D. 


[Corollary 3] The optimal prediction estimate u( t,x/ t^ ) and prediction error 
covariance function P(t,x,y/t. ) can be represented as 


u(t,x/tj = / GU,t 1/ ,x,a)u(t k ,a/t |< )da 


^ / J ^ ! 


(27) 


and 


P(t,x,y/tJ = // 6( t , tr , x ,ot) P ( t^ ,cx , 3/ 1^ )G( t , t^ ,y , B) 
k D D K 


dadp 


*fu. G(t,o,x,a)Q(a,a,B)G(t,a,y,B) dadBda. 


(28) 


t k DD 


[Proof] It is clear that (19) and (22) possess unique solutions. Differen- 
tiating (27) and (28) with respect to t yields (19) and (22), 
respectively. Since (19) and (22) have unique solutions, (27) and (28) are 
those solutions. Q.E.D. 
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V. Derivation of the Optimal Filter 

In order to derive the optimal filter by using the Wiener-Hopf theorem 
for the filtering problem, we represent the solution of (3) in terms of the 
fundamental solution G(t,a,x,y) as 


u(t k+i ,x) = / G (t k+1 »t k ,x,a)u(t k ,a) da 


f ^ +1 f 

+ J J G(t + 1 ,n,x,a)w(ri,a) 


t k D 


k+r 


dadn 


(29) 


and 


u t k+ 1 ^ x 3 ^ " l VVl’VV®^ V 0 ^ da 


r k+1 r 

+ J J G M (t k+ 1 ,n>x 3 a)w(n,a) dadn 


*k ° 


(30) 


where 


G M^ ^k+'i ,n ’ x 3 ,ct ^ 


s{t k+1 .n, x 1 1 (k t 1) ,x^ k+1) ,x 3 ,o) ' 
G(t kn ,n, xf ;i ),xf +1 ),x 3> a) 


(31) 


From (17) we have 


F ( t k+ 1 ,x,t k+ 1 )EtZ(t k+l) Z " (t ?)' 1 


+ 22 F(t k+1 ,x,t 0 )ElZ(tJZ"(t r )] = E[u(t m ,x)Z^(t r )] 


a=0 

for ; = 0,1 , , k+1. 


k+1’ 


(32) 
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From (29) and the independence of Z(t^), t, - 0,1,..., k and w(ri,x),. 
t k < n < t^ +1 it follows that 


E[u(t k+1 ,x)Z'(t r )] = / G(t u+1 ,t^,x,a)E[u(t u ,a)Z'(t r )] 


k+1 k : 


da. 


Using the Wiener-Hopf equation (17), we have 


Etu(t k+1 ,x)Z'(t )] 


/ G(t 


k + i,t k ,x,a) 


E 

c=0 


F(t k ,a,t a )E[Z(t a )Z"(t ? )3. 


On the other hand, from (8) and the whiteness of v(t k+1 ), we have, for 

^ V 


E[Z(t k+ i)Z'(t c )] 


rh 

/ J(t 


k+1 ,x 3 )E[u t 


k+1 


(x 3 )Z"(t )] dx 3 . 


Substituting (30) into the above equation and using the independence of 
Z(t ), t < t k and w(n,a) , t R < n < t k+J yields 

E[Z(t k+1 )Z"(t c )] = / h J(t k+1 ,x 3 ) _[ G M (t k+1 ,t k ,x 3 ,a)E[u(t k ,a)Z"(t c )] dadx 3 

Again, we use the Wiener-Hopf equation (17) in the above equation and 

h k 

E[Z(t k+1 )Z'(t ? )] = / J(t k+1 .x 3 ) | G H (t k+r t k ,x 3 ,c,)'i:F(t k .a.t 0 ) 

E[Z(t a )Z"(t )J dadx 3 . (34) 

Substituting (33) and (34) in (32) yields 
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E f 4<V Xlt o )E - 0 

where 

F A (t k ,x,t 0 ) = F ( t | <+ i» x » t | < +i) f J ( t k+r x 3^ G M^ t k+r t k 5X 3’ a ^ F ^k’ a ’ t o^ dadx 3 
+ F(t k ^>x,t^) " j G( j x s ct ) F( ,cx s t^) dot. 

Then from Lemma 1 we have F^{t^,x,t c ) = 0, and we have the following lemma. 
[Lemma 2] The optimal kernel function F(t k+1 ,x,t a ) of the filter is given by 

F(tk+l ’X»t 0 ) - j 8 ( ^k+l ’ ^k * x F ^ ^k d ^ 

- F(t k+ i ,x 9 t k+1 ) j J * t k +i’ x 3 ) / G M (t k+l ,t k ,x 3’ a)F(t k ,a,t o ) dadx 3' 

[Theorem 4] The optimal filtering estimate u (^k+l’ x ^k+l^ 1S 9 lven 
0(t k+r x /t k+ i) ■ 0(t k+1 ,x/t k ) + F ( t k+i’ Xit k+l^ v ^ t k+l^ 

v(t k+ i) k Z(t k+1 ) - / d(t k+1 ,x 3 )u t ^^(x 3 /t k ) dx 3 , (37) 

u(t a ,x/t 0 ) = u Q (x), (38) 

r ^(tk +1 ,£/tk +1 ) = 5(t k+1 ,5), C £ 9D (39) 


where 
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u t (x 3 /t k ) = 

\+l J K 


G <Vr x } (k : 1) ' x 2 (k+1> - x 3 /t k ) 

u(t x M(k+1 > x M(k+1 > X /t ) 

y lt k+l’ x l ,x 2 * X 3'V 


[Proof] Using (12) and (35) yields 


“ (t k+r x /W = F ^ t k+r x,t k+P z ^ t k+P 
, K 

+ D 6(t k+l» t k* x * ft) F( V a,t 0 )Z(t a> da 


h ■ 

- F ( t k+ I> x » t | <+ 1 ) / J(t k+l’ x 3 ) I G M^ ^k+l * t k * x 3 * a ^ ^ F (t k >a«t )Z(t ) dadx 3' 
u v 0=0 


Then from (12) end (27) we have 


^ j * x/ t^i ) f S ( j $ X > ex) u ( > ot/^^ ) da 
r h r 

+ F ( t k +l ,x,t k+l^ Z ^ t k+l^ ’ J J(t k+l ,x 3^ J Q G M^ t k+l * t k ’ x 3 ,a)u(t k> a/t k ) dadXg 

= u(t k+1 ,x/t k ) + F (t k+1> X,t k+1 )v(t k+1 ). 

Since the initial and boundary conditions are clear, the proof of the theorem 
is complete. Q.E.D. 

To determine the optimal kernel function F (t k+1 >x,t k+1 ) , we introduce 
the following notation, 

P M (t x ,x,y 3 /t k ) = [P(t T ,x,y 1(k) /t k ) P(t T ,x,y M(k) /t k )] (40) 


and 
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WV x a’W 


P(t T ,xJ {k) s y 3 /t k ) 

* 

P(t T ,x M(k) ,y 3 /t k ) 


(41) 


P(t T ,x^ k \y^ k Vt k ), .... P(t T ,x^ k ^y M ^ k Vt k ) 

.... P(t T ,x M ( k ),y M ( k )/t kl 

where x m < k > = (x* (k) , x" 2 l(k) , x 3 ) and y m ( k > = (yf k \ yf k K y 3 ). m - 1,2,..., M. 
From the definitions of P^(t T ,x,y 3 /t k ) and P^ M (t T ,X 3 ,y 3 /t k ) it follows that 


and 


P M (t T ,x,y 3 /t k ) = E[u(t T ,x/t k )Q t (y 3 /t k )l 




(42) 


(43) 


where 


and 


Q t (x 3 /t k ) = u t (x 3> ‘ G t ( W 


u tJW s 


u(t T ,x^ k Vt k ) 
G(t T ,x M(k) /t k ) 


(44) 


(45) 


Furthermore, we define the covariance matrix of the innovation process v (t k+ j) 
by r(t k+ 1 /t k ) = E[v(t k+1 )v'(t k+1 )]. Then from (37) we have 


-h rh 


r ^k+l^V S I / J ^ t k+1 ,X 3^ P MM^ t k+l * x 3 ,y 3 /,t k^ J ^k+l’ y 3^ dx 3 dy 3 


o o 


+ «W 


(46) 
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[Theorem 5] The optimal filtering gain function F ( is 
given by 


r h , 

F(t k+1 ,x,t k+ i) J P M ( t k+r x ,y 3 / t k ) 0 " ( t k+l ,y 3 ) dy 3 r (t k+1 /t k ). (47) 

[Proof] From the Wiener-Hopf equation (17) we have 

F (t k+ i,x,t k+1 )E[Z(t k+ i)Z"( t k+ i)] 

k 

+ E q = E[u<t k+1 ,x)Z'(t k+1 )]. 

Substituting (35) into the above equation yields 

F ( t k+l * x » t k+l ) E t (Z( t k +i) " / J ( t k+l’ x 3^t k+1 ^ x 3 /t k^ dx 3^^ t k+l^ 

= F [ (u( t k _j_^ ,x) - u(t k _j_j >x/t k ) )Z (t^j)). 


Using (8) and the orthogonality condition of Corollary 1 yields 


m\+ v x/t k )Z'{t k+l )] ~ f E[u(t k+ i,x/t k ) u J (x 3 )]J (t k+1 » x 3 ) dx 3 

o k+i 

h 

P M^ t k+l ,x ’ x 3 /t k^ J "^ t k+r x 3^ dx 3 



and 



,x 3^ P MM^ t k+r X 3 ,y 3 //t k^ J ^ t k+l’ y 3^ dx 3 dy 3 


+ R( W = r <WV- 


(48) 
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Then we have 

F(t k+l ,x,t k+l )r(t k+l /t k ) = f P M ( Vl> x ’ x 3 /t k ) ' r(t k+l-’‘3 ) dx 3 
and the proof of the theorem is complete, Q.E.D. 

[Theorem 6] The optimal filtering error covariance function P(t k+1 ,x,y/t k+1 ) 
is given by 


P(t k+1 ,x, y/ t k+ i) = P(t k+1 ,x.y/t k ) 


rh rh 


/" /" P M (t k+1 ,x,X3/t k )J1t k+1 ,x 3 )r- 1 (t k+! /t k )J(t k+1> y 3 )P | J(t k+1 ,y,y3/t 1( ) dxjdyj 


0 o 


P(t 0 ,x,y/t 0 ) = P 0 (x,yj 
r e P(t k+1 ^y/t k ) = o, £ € 3D. 

[Proof] From (3) and (36) we have 

G(t k+1 ,x/t k+1 ) = u(t k+1 ,x/t k ) - F (t k+1 »x,t k+1 )v(t k+1 ) 

and from (4) and (39), 


(49) 

(50) 

(51) 


(52) 


r e u(t k+1 ,5/t k+1 ) - o. 


£ € 3D. (53) 

Using the independence of v(t k+1 ) and u(t k+1> x/t k ) or ^k+l’^k^ 
P(t k+1 ,*,y/t k+1 ) « E[«t k+r x/t k+ pn(t k+1 ,y/t k+1 )] 
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= P(t k+1 ,x,y/t k ) + F(t k+1 ,x,t k+1 )E[v(t k+1 )v"(t k+1 )]F"(t k+1 ,y,, v . k+1 ) 

" ^k+l ,x,t k+l^ / ^k+1 ’ x 3^^t u , , ^S^k^^k+l’^^k^ ^ x 3 
0 K+l 

- / E[Q(t k+1 ,x/t k )Qj (y 3 /t v ))r(\ +v y 3 )dy 3 r(t k+v y,t m ). 
o k+1 

Using (40) and (47) yields 

P(t M ,x,y/t k+1 ) - P(t k+1 .*.y/t k ) 

h h 

" / / *V t k+l’ x ’ x 3 //t k^ J (t k+l ,x 3 )r ^k+l^k^ 
o o 

J(t k+l’ y 3 )P M (t k+l’ y ’ y 3 /t k ) dX 3 dy 3- 

Since the initial value u(t Q ,x/t 0 ) is equal to u o (x), it is clear that 
P(t 0 »x,y/t 0 ) = P Q (x,y). Multiplying each side of (53) by Q(t k+1 ,y/t k+ i) and 
taking the expectation yields (51). Q.E.D. 
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VI. Derivation of the Optimal Smoothing Estimator 

In this section we derive the optimal smoothing estimator by using the 
Wiener-Hopf theory. 

[Lemma 3] The optimal kernel function B(t T ,t k+ j,x,t 0 ) of the smoothing estima- 
tor is given by 


B(t T ,tk + i >x,t a ) = B(t^ . t^ »x , t 0 ) 

- B(t T ,t k+r x,t k+1 ) / J(t k+1 ,x 3 ) / G M (t k+1 ,t k ,x 3 ,a)F(t k ,a,t 0 ) dadx 3 . (54) 


[Proof] From the Wiener-Hopf equation (18) for the smoothing problem we have 


k+1 

^ B(t Tj t k+1 ,x,t u )E[Z(t o )Z'(t ? )] = E[u(t T ,x)Z"(t c )], 
5 = 0,1,. .. , k+1 


(55) 


and 


k 

^ B(t T ,t k ,x,t o )E[Z(t 0 )Z'(t c )] = E[u(t T ,x)Z"(t ; )], 
K = 0,1,..., k. 


(56) 


Subtracting (56) from (55) yields 


B ^ t x ,t k+r x ’ t k+l^ E ^ Z ^ 1 : k+l^ Z 


k 

+ £ (B(t T ,t k+1 .x.t 0 ) - B(t T ,t k ,x.t 0 ))E[Z(t a )Z'(t )1 = 0. 
o=0 


From (8) and (17) we have 
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F[Z(t k ^)Z (tj-)3 = J ’^3^f 0 ^ ( ^k+l * 3 °t) E [u ( tj^ 9 cx) Z (t^,)] dcxdx^ 

^ o 0 

h - k 

= / J ^k+r X 3^i G M^’ t k+l*' t k* x 3* <x ^ 5 V Z ^V ] dadx 3* 

0 U 0~U 

Then it follows that 
k 

X F A (t T .t k ,x,t o )E[Z(t a )7'(t )] ■ 0 
a = 0 

where 

F/ ( tj J t|^ 9 X 9 t Q ) = B ( ^ 3 ^q) " B( t^) 

+ / J (Vl’ x 3 ) I G M^ t k+r t k ,x 3 ,a ^ tr ^ t k’ a,t a^’ 

0 u 

Since it is clear that B(t ,t k ,x,t 0 ) + F A (t T> t k ,x,t a ) also satisfies the Wiener 
Hopf equation (18), from Lemma 1 F A (t , tk» x »t 0 ) = 0, a = 0,1,..., k. Thus, 
the proof of the lemma is complete. Q.E.D. 

[Theorem 7] The optimal smoothing estimate uCt^x/t^) is given by 

u(t T ,x/t k+1 ) = u(t T ,x/t k ) + B ( V Vl’^k+l^^k+l^ ( 57 ) 

= s(Tj ^ 5 5 £ 9D - (58) 

[Proof] From (13) it follows that 

u(t T , X /t k+i ) = B( » x ,t k+1 )Z( t k+1 ) 
k 

+ ^VVl’^V^ V’ 
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Substituting (54) into the above equation yields 


u(t T> x/t k+1 ) = B ( V^+l’^k+l^^k+l^ 


+ 1 B(t ,t k> x,t 0 )Z(t 0 ) 
o=0 


and substituting (13) into the above equation yields (57). Since we have no 
additional information about the boundary value of u(t T ,x) except for S(t x ,f;), 
we have (58). Thus, the proof of the theorem is complete. Q.E.D. 

[Theorem 8] The optimal smoothing gain function B(t T ,t k+1 ,x,t^ +1 ) is given by 

r h -i 

B( Wl’ X ’W = J N ( t x » x » x 3/ t k+p J ^k+r X 3^ dx 3 r ^k+l^k^ 

o 

where 


N(t T ,x,Xg/t k+ j) = / M(t^,x,y/t k )G^(t| <+ 1 ,t k ,x 3 »y) dy 


M(t x ,x,y/t k ) = E[u(t T ,x/t k )u(t k ,y/t k )]. 

[Proof] From the Wiener-Hopf equation (18) we have 

B ^ t T ,t k+r x,t k+l^ E[Z ^ t k+l^ Z ^k+l^ 

k 

+ 1 B(t T ,t k+1 ,x,t o )E[Z(t a )Z^(t k+1 )] = E[u(t T ,x)Z"(t k+1 )]. 
g=0 


(60) 


(61) 
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Substituting (54) into the above equation yields 


^ ( t*c * tk+i 1 x ’ tj <+ i ) E [\>( Lk + i )Z (t k+ ,)] “ E [u( t T ,x/t^) Z U k+1 )] 


On the other hand, from (27) and (29) 


( 62 ) 


Q(t k+ i,x/t k ) = / G(t k+1 ,t k ,x,y)u(t k ,y/t k ) dy 

A+lf 

+ J J G(t k+1 ,n,x,y)w(n,y) dydq. 


t k D 


Then we have 


E[u(t x ,x/t k )Z'(t k+1 )] = / L M ^ t T ,x,y / t k) G M( t k +r t k’ x 3 ,y ) dy J ^k+l ,x 3 ) dx. 


o D 


r h 

J N(t_,x,x 3 /t k+1 )J ( 9 x 3 ^ dx : 


Substituting (48) and the above equation into (62) yields (59). Thus, the 
proof of the theorem is complete. Q.E.D. 

Let us now derive the equation for M(t x ,x,y/t k+1 ) . Using the orthogonality 
condition of Corollary 1 yields 


M(t T ,x,y/t k+1 ) = E[u(t T ,x)u(t k+1 ,y/t k+1 )] 
Substituting (52) into the above equation yields 


(63) 


M(t T ,x,y/t k+1 ) = / G(t k+1 ,t k ,y,a)M(t T ,x,a/t k ) da 


rh 

- J N(t T ,x,x 3 /t k+1 )iT(t k+1 ,x 3 ) dx 3 F 
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r g a(t T ,?/t k+1 ) = o, 5 € 3D. 


( 64 ) 


Multiplying each side of the above equation by u(t 
expectation yields r^MU^^y/t^) = 0, 5 € 3D. 
holds. 


k+i »y/ t k+i ) and takin s the 
Thus, the following theorem 


[Theorem 9] M(t T ,x,y/t k+1 ) is given by 


M(t^,x,y/t k+ j) - j G (t k+1 >t k ,y,a)M(t T ,x,a/t k , do 


f h 

J N(t T> x,x 3 /t k+1 )J"(t k+1 ,x 3 ) dx 3 t k+1 » 

(65) 

M(t T ,x,y/t, r ) = P(t T ,x,y/t T ),. 

(66) 

r ? M(t T ,?,y/t k+1 ) = 0, 5£3D. 

(67) 


It remains to derive the equation for the optimal smoothing error covariance 
function P(t ,x,y/t k+1 ) . From (57) we have 

u(t T ,x/t k+1 ) = u(t T ,x/t k ) - B(t T ,t k+1 ,x,t k+1 )v(t k+1 ). (68) 

[Theorem 10] The optimal smoothing error covariance function P(t »x,y/t k+ ^) 
is given by 

P(t T ,x,y/t k+1 ) = P(t T ,x,y/t k ) 

"If N ( t x » x » x 3/ t | c+ i^''f t k+r x 3^ r ^k+l^k^Vt-l'^ 

0 0 

N( t T .y.y 3 /t k+1 ) dx 3 dy 3 , 


(69) 
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r ? P(t T ,e,y/t k+1 ) =0, 5 € 3D. (70) 

[Proof] From ( 68 ) we have 

P(t T ,x,y/t k+1 ) = E[D(t T ,x/t k+ 1 )u(t T ,y/t k+1 )] 

P(t T ,x,y/t k ) + B ( t T » t | <+ 1 >x>t 1 <+ 1 )r(t k+ 1 /t k )B"(t T ,t 1 <+ 15 y,t l<+1 ) 

- B ( t T ’ t k+l ,x ’ t k+l )Etv,(t k+l )a(t T’ y/t k )] 

- E[0(t T ,x/t k )v"(t k+ 1 )]B"(t T ,t k+ 1 ,y,t k+1 ). (71) 

But we have 

E[Q(t T ,x/t k )v'(t k+1 )] = // G^(t k+ 1 ,t k ,x 3 ,a)M(t T ,x,a/t k ) 

o D 

0"(t k+ i * x 3 ) dadx 3 


E[v( t k+ i)u(t T ,x/t k )] 


f h ! 


J ( t k _^ 5 ^ 3 ^ B |_ ( ^k +1 ?t k ,X 3 ’ a ^ ^ dadx 


3* 


Substituting the above equations and (47) into (71) yields (69). Multiplying 
each side of (64) by u(t T> y/t. k+ ^) and taking the expectation yields (70). Q.E.D. 

[Theorem 11] The optimal smoothing estimator is given by 

k 

u(t T ,x/t k ) = u(t x ,x/t T ) + J B(t T ,t Jl ,x,t £ )v(t Jl ) 


(72) 
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and the optimal smoothing error covariance function P(t T ,x,y/t k ) is given by 
P(t T ,x,y/t k ) = P(t T ,x,y/t T ) 

' J -,+1 / h / h|,(t f* ,X 3 /t * ) ' , ' (t * ,X 3 )r ‘ 1( V t *-l ) 

J (t Jt .y 3 )N(t T .y,y 3 /t i ) dx 3 dy 3> (73) 
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VII. Estimation of the Concentration Distribution Downwind of a Continuous, 
Ground-Level Line Source 

There has been much recent interest in the airborne measurement of pollu- 
tant concentrations downwind of sources [9 ] - [11], Here we wish to consider 
a hypothetical, but realistic, situation in which an aircraft with a downward- 
looking instrument, such as for example the JPL Laser Absorption Spectrometer 
[12], is flown at different altitudes downwind of the source, and total species 
burdens are measured at a series of downwind distances. 

The steady-state concentration of a species downwind of a continuously 
emitted ground-level line source (e.g. a highway) situated normal to the 
direction of the wind flow is governed by the following form of the atmospheric 
diffusion equation [5], 


V 


3u_ 

i ax. 


JL 

ax’ 


k v (x 3 ) 


'i 3 \ 
u (0,x 3 ) = u 0 (x 3 ) 


K v (0) = *«( Xl ). 



au 

9x 3 


+ w(xj,x 3 ) 


x 3 = 0 
x 3 = h 


(74) 

(75) 

(76) 

(77) 


where <j> is the constant rate of release. For convenience we will take K y = 1, 
since vertical variations of this constant are not essential to the estimation 
problem we will consider. If we let t = x^/V^ and x = x^, (74)- (77) become 


I = 4 + «<*.*> 

oX 

u(0,x) = Uq(x) 


(78) 


(79) 
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= - 4>6(t), x = 0 


x = h. 


In this case the measurements ZU^) are related to the concentration 
u(t k ,x) by (8), 

r h 

Z(t.) = J J n (t.,x)u(t ,x) dx + v ( t k ) 


where the instrument kernel function will be taken to have the form, 


0"(t k ,x) * 




x > h. 


n 1,2, , 


The theory developed in the prior sections can be applied directly to this 
problem, and the optimal filter and smoother are given in Table 1. The pre- 
diction, filtering and smoothing algorithms were applied to hypothetical data 
generated by solving ( 74 ) - (77 ) and forming Z(t^) from (82), using noise 
processes w(t,x) and v( t^) with prescribed properties. The algorithms were 
applied to estimate the concentration distribution u(t k ,x) as a function of 
height x at several downwind distances, t^, t£,... based on measurements taken 
at one to four elevations. It is of interest to study the behavior of the 
estimates as a function of downwind distance and of the number of elevations 
at which data are simultaneously taken. Values of all parameters used in the 
calculation are given in Table 2. 
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Figs. 2-4 show selected results of the application of the filtering and 
smoothing algorithms to the synthetic data of this example. Fig. 2 shows a 
comparison of the true concentration distribution u(t 1} x) and the filter esti- 
mates, G(tj,x/tj) based on two and four measurement elevations (tj = 0.0002). 

As expected, the profile estimated on the basis of four measurement elevations 
is superior to that based only on two altitudes. Fig. 3 shows similar results 
at tg s 0.0082. The filter estimate based on n = 4 virtually coincides with 
the actual concentration distribution. The performance of the smoothing algo- 
rithm is illustrated in Fig. 4, in which the true concentration u(t ,x) is com- 
pared with the filter estimate, u(t T ,x/t T ), and the smoothed estimates, 
u(t T ,x/t 2 ), and u(t T ,x/t^), with t T = 0.0002, t 2 = 0.0012, and t^ = 0.0032. 

Table 3 gives the trace of the filtering error covariance matrix, P(t,x,x/t), 
for the four measurement configurations at three downwind distances t. As expec- 
ted, the trace decreases as the number of measurement elevations is increased 
from 1 to 4. 

VIII. Conclusions 

Filtering and smoothing algorithms for the processing of remote sensing 
data on atmospheric species concentrations have been derived using Wiener-Hopf 
theory. The algorithms were applied successfully to estimate concentration 
distributions from a hypothetical ground-level line source of material (e.g. a 
highway) based on remote sensing data taken from several elevations at a number 
of points downwind from the source. Although there has been increasing interest 
in the remote sensing of airborne concentrations, a data set sufficient for ap- 
plication of the theory developed in this paper does not yet appear to exist. 
Nevertheless, it is hoped that the availability of the algorithms developed here 
will facilitate processing of remote sensing data in conjunction with mathematical 
models of air pollutant behavior. 
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Table 2. Parameter Values Used in Line Source Estimation 

Example 


Truncation number N = 5 


Measurement time t k+J = t k + 0.001, k = 1,2,3,*** 
where t-j = 0.0002 

Fixed-point time for smoothing t = 0.0002 
Constant rate of release 4> = 0.3 

Measurement points h p = n/4, n = 1,2,3, 4 
Initial values and noise covariances 

Etu (x)] = u (x) = J ufo(x) 

U u i=1 1 1 

N n 

Covfu 0 (x),u 0 (y)] = P Q (x,y) = J (x)<j>. (y) 


Cov[w(t,x) ,w(s,y)] = Q(t,x,y)6(t-s), Cov[v(t k ),v(t n )] = R( t k )6 kn 
N 

Q(t,x,y) = .I q 1 ^ 1 (x)^(y), R(t k ) = d1ag[r^ ,r 2 ,r 3 ,r 4 ] 


^ (x) 


1 

/2 cos(i-l )ttx 


X, - - (i-l)V 


i = 1 

1 > 2 
i > 1 


i 

i 

2 

3 

4 

5 

“? 

3.0 

1.0 

0.03 

0.003 

0.0003 

p?i 

< 

0.1 2 

0.01 2 

0.001 2 

0.0001 2 

q ii 

1 

0.5 2 

0.25 2 

0.125 2 

0.0625 2 

r i 

0.1 2 

0.07 2 

i 

0.05 2 

0.03 ? 

-——————-pr 
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Table 3. Trace of the Filtering Error Covariance Matrix P(t,x,x/t) 


Measurements 

t = 0.0002 

t = 0.0032 

t = 0.0062 

4 point 
( h-j ^4) 

0.2405X10" 1 

0.11 89x1 0 -1 

0.9616xl0' 2 

3 point 
( h -1 5 h ^ > h^ ) 

0.3441 xl 0"^ 

0 . 1577x10 

0,1 335x1 0" 1 

2 point 
(h r h 2 ) 

0.1678 

0.1195 

0.1137 

1 point 
(h } ) 

0.6700 

0.2914 

0.2338 
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Figure Captions 

Fig. 1. Remote Sensing Measurement Configuration Considered in This 
Work. 

Fig. 2. Comparison of True Concentration u(t|,x) and the Filter Esti- 
mates u(t|,x/t^) based on 2 and 4 Measurement Elevations. 
t ] = 0.0002. 

Fig. 3. Comparison of True Concentration u(tg,x) and the Filter 

Estimates u(tg,x/tg) based on 2 and 4 Measurement Elevations. 

tg = 0.0082. 

Fig. 4. Comparison of True Concentration u(t ,x), the Filter Estimate 
u(t ,x/t ) and the Fixed Point Smoothing Estimates u(t T ,x/t 2 ) 
and u(t ,x/t 4 ). t = 0.0002, t 2 = 0.0012, t 4 = 0.0032. 
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Figure 1 
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Figure 4 



SUMMARY AND CONCLUSIONS 

The object of this research grant was to initiate an evaluation of the 

analysis of remote sensing data on pollutant concentrations in the troposphere. 

Remote sensing measurements of pollutant concentrations are becoming increas" 

ingly important in understanding the transport and transformation of pollutants 
* 

over moderate to long distances in the atmosphere. Traditionally such data 
have not been analyzed beyond the point of constructing mass fluxes and total 
budgets over a region. The question studied in this research grant was that 
of the further analysis of such data, particularly when one has a mathematical 
model available. The specific problem then is to see how typical remote sens- 
ing data can be used in conjunction with a mathematical model to extract addi- 
tional information about the pollutant behavior in the region being studied. 

The essential problem is one of estimation, that is, of using the typical 
remote sensing data to determine full concentration distributions. Once full 
concentration distributions are available, one can then assess the mechanisms 
of the process through the mathematical model. The first step in the research 
was to look theoretically at the question of the minimum amount of data needed 
to reconstruct a concentration distribution from finite data typical of those 
collected in remote sensing. Chapter I of this report presents a development 
and derivation of a condition of reconstructability, namely rigorous conditions 
that can be applied to a data sampling program to determine whether it will be 
possible to estimate a species concentration distribution from such measure- 
ments. Chapters II and III of this report are then devoted to the development 
of a numerical algorithm that will process the data to produce concentration 
distribution estimates in the cases when the data are a priori reconstructable. 

Perhaps the most important result of this study is the indication of the 
types of measurement strategies one should follow in remote sensing programs. 


In particular, it appears that the best measurement strategy is to attempt to 
obtain pollutant burdens at a certain location at a number or elevations at 
times as close as possible. This strategy is recommended because the vertical 
distribution of pollutant concentrations in the first 1,000 meters of the atmo- 
sphere is a crucial element of a mathematical model of such species. The 
theory and numerical techniques developed in this study will tell one when 
devising a measurement program and monitoring strategy the number of vertical 
levels at which one should make measurements to be able to estimate relatively 
accurately the complete vertical concentration profile of the species of inter- 
est. It is anticipated that these results will be of value to those contem- 
plating remote sensing measurement programs of tropospheric species that involve 
measurements at several vertical levels. 


